Mendelevium

Diary
Drug Design
Field Knowledge
- Academia
  - Yang
- Biology
- Physics
Free Energy
Machine Learning & AI
Molecular Dynamics
Other
Specific Sytems
Techniques
- Linux
- Python
- Research
- Web
about
Home

Contact

Copyright © 2025 Xufan Gao | Academic Research Blog

Home > Molecular Dynamics > Sampling & Analysis > Allostery

A Bunch of Biophysics is Loading ...

Allostery

神经关系推断：从MD轨迹中学习蛋白质长程变构相互作用
 神经关系推断：从MD轨迹中学习蛋白质长程变构相互作用本文信息标题：Neural Relational Inference to Learn Long-range Allosteric Interactions in Proteins from Molecular Dynamics Simulations 作者：Jingxuan Zhu¹,²,³, Juexin Wang¹,², Weiwei Han¹, Dong Xu² 发表时间: 2022年3月10日单位：吉林大学生命科学学院，酶学与工程教育部重点实验室（中国长春）密苏里大学电气工程与计算机科学系，Bond生命科学中心（美国哥伦比亚）期刊：Nature Communications 引用格式：Zhu, J., Wang, J., Han, W. & Xu, D. Neural relational inference to learn long-range allosteric interactions in proteins from molecular dynamics simulations. Nat Commun 13, 1661 (2022). https://doi.org/10.1038/s41467-022-29331-3 源代码：https://github.com/juexinwang/NRI-MD 摘要蛋白质变构是一种由空间上长程的分子内通信促进的生物过程，即远端位点的配体结合或氨基酸变化能够远程影响活性位点。分子动力学（MD）模拟为探测变构效应提供了强大的计算方法。然而，当前的MD模拟仍无法达到整个变构过程的时间尺度。深度学习的出现使评估空间上短程和长程通信以理解变构成为可能。为此，我们应用了一种基于图神经网络的神经关系推断模型，该模型采用编码器-解码器架构同时推断潜在相互作用，将蛋白质变构过程探测为相互作用残基的动态网络。从MD轨迹中，该模型成功学习了可以介导Pin1、SOD1和MEK1系统中远端位点间变构通信的长程相互作用和路径。此外，该模型能够在MD模拟轨迹中更早发现与变构相关的相互作用，并比其他方法更准确地预测突变后的相对自由能变化。核心结论深度学习破解变构难题：首次将神经关系推断（NRI）模型应用于MD数据分析，通过encoder-decoder架构从MD轨迹中推断残基间的相互作用网络长程通信路径识别：成功识别了Pin1、SOD1和MEK1三个系统中介导变构通信的长程路径，揭示了WW域与催化位点之间的通信机制早期信号捕获能力：NRI模型能在MD轨迹的早期阶段（50-100 ns）检测到变构信号，远早于传统方法（200 ns以后）自由能预测优势：基于学习到的相互作用网络计算的自由能变化与实验数据高度一致（$R^2=0.939$），显著优于传统方法（$R^2=0.188$）物理可解释性：学习到的相互作用类型具有明确的物理意义，揭示了结构域间的动态耦合模式背景蛋白质变构是蛋白质功能调控的核心机制之一，通过空间上远离活性位点的区域（如别构位点）来影响蛋白质的活性。这种长程通信机制使蛋白质能够整合多个信号输入，实现精细的功能调控。然而，理解变构信号如何在蛋白质内部传播一直是结构生物学领域的重大挑战。传统研究变构的方法主要基于静态晶体结构或简化的弹性网络模型，但这些方法难以捕捉蛋白质在全原子模拟中的动态复杂性。分子动力学（MD）模拟虽然能够提供原子级别的运动信息，但由于变构过程通常发生在微秒到毫秒时间尺度，而常规MD模拟仅能达到纳秒到微秒级别，使得直接观测完整的变构过程变得困难。近年来，图神经网络（GNN）在分析复杂系统方面展现出巨大潜力。特别是神经关系推断（NRI）模型，作为一种无监督学习方法，能够同时推断系统中实体间的相互作用关系并预测系统演化。这种方法已被成功应用于交通系统、动态物理系统和计算机视觉等领域，但在生物分子系统中的应用尚属空白。关键科学问题时间尺度不匹配：MD模拟的时间尺度（纳秒-微秒）远短于完整变构过程（微秒-毫秒），如何从有限长度的轨迹中提取有意义的变构信息高维数据分析困难：MD轨迹产生的高维（$3N$维）动态数据难以直接分析，需要有效的降维和信息提取方法因果vs相关关系：传统基于相关性的方法难以区分变构通信中的因果关系，可能误判非因果性的相关关系长程通信识别：如何在复杂的残基相互作用网络中准确识别介导长程变构通信的关键路径创新点 NRI模型首次应用于MD分析：首次将神经关系推断模型应用于生物分子MD数据分析，通过GNN同时推断残基间的潜在相互作用动态相互作用网络：将蛋白质变构过程建模为相互作用残基的动态网络，学习到的边权重反映了残基间相互作用的强度轨迹重建验证：通过重建原始MD轨迹来验证学习到的相互作用的有效性，确保模型捕获的是真实的物理相互作用早期信号检测：NRI模型能够在MD轨迹的早期阶段（50-100 ns）检测到变构信号，比传统方法提前数倍自由能准确预测：基于学习到的相互作用网络计算突变后的相对自由能变化，与实验数据高度一致研究内容 NRI模型架构与训练图1：通过重建MD模拟轨迹推断相互作用图的过程该图展示了NRI模型的完整工作流程，从系统准备到相互作用推断： (a) 变构系统准备：准备配体-结合复合物或突变蛋白质的变构系统结构，包括Pin1（WW域+PPIase域）、SOD1（β桶+活性环）、MEK1（N叶+C叶+激活片段） (b) MD模拟：对制备的变构系统进行MD模拟，获得包含动态3D坐标的轨迹数据，采样间隔约为20 ns，总模拟时间100-500 ns (c) 常规分析：传统的MD轨迹分析方法，如RMSD、RMSF、PCA等，提供结构变化和柔性信息 (d) NRI模型：包含两个 jointly 训练的组件——编码器（推断潜在相互作用的因子化分布$q_\phi(z x)$）和解码器（基于采样的相互作用重建动态系统）编码器-解码器架构 NRI模型的核心思想是将MD轨迹中的残基运动建模为动态系统，其中每个残基的运动受到其与其他残基相互作用的影响。模型采用变分自编码器（VAE）框架，最大化证据下界（ELBO）： \[\log p_\theta(x) \geq \mathbb{E}_{q_\phi(z|x)}[\log p_\theta(x|z)] - D_{KL}(q_\phi(z|x) || p_\theta(z))\] 其中： $x$ 是MD轨迹中的残基坐标 $z$ 是残基间的潜在相互作用（以边的形式表示） $q_\phi(z x)$ 是编码器推断的后验分布 $p_\theta(x z)$ 是解码器重建的轨迹分布 $p_\theta(z)$ 是先验分布（均匀独立的分类分布）编码器采用图神经网络（GNN）在完全连接网络上处理输入坐标，输出每个残基对的相互作用类型分布： \[q_\phi(z_{ij}|x) = \text{softmax}(f_{\text{enc},\phi}(x)_{ij,1:K})\] 其中 $K$ 是相互作用类型的数量（本文中$K=10$），$f_{\text{enc},\phi}(x)$ 是GNN编码器。解码器根据采样的相互作用$z$重建动态系统，预测下一时刻的残基位置。通过最小化重建误差（MSE）和最大化似然，模型学习到有意义的相互作用模式。 GNN消息传递机制：Receive与Send NRI模型的核心是图神经网络的消息传递机制，通过交替的”节点到边”和”边到节点”操作来传播信息：节点到边（Send）操作：节点发送自身嵌入给相连的边对于每条边$(i,j)$，接收来自节点$i$和节点$j$的嵌入： \[h_{ij} = f_e([h_i, h_j])\] 物理意义：节点向可能的相互作用伙伴传达自身状态信息，这里$h_i$和$h_j$是节点的隐藏状态表示。边到节点（Receive）操作：节点接收来自所有连接的边的消息节点$j$接收的消息： \[h_j^{\text{new}} = f_v\left(\sum_{i \neq j} h_{ij}\right)\] 物理意义：节点整合来自所有相互作用伙伴的信息，更新自身的状态表示。这里$\sum_{i \neq j} h_{ij}$表示聚合所有指向节点$j$的边消息。多轮消息传递：初始节点嵌入：将轨迹特征映射到节点嵌入$h_i = f_{\text{enc}}(x_i)$ 第一轮v→e：计算所有残基对的边嵌入候选$h_{ij}$ 第一轮e→v：聚合边消息更新节点状态重复：进行多轮消息传递（通常2-3轮）生成分布：输出每条边的$K$种相互作用类型分布$z_{ij}$ 这种机制使模型能够捕获残基间复杂的、非线性的相互作用模式，而非简单的线性相关或距离依赖关系。 graph TB Start["MD轨迹输入<br/>N个残基×T帧×3维坐标"] --> Encoder["编码器 (GNN)<br/>推断相互作用z_ij"] Encoder --> Latent["潜在变量<br/>z_ij ∈ {1,...,K}<br/>K种相互作用类型"] Latent --> Decoder["解码器<br/>重建轨迹x'"] Decoder --> Loss1["重建损失<br/>MSE(x, x')"] Encoder --> Loss2["KL散度<br/>正则化先验"] Loss1 --> Joint["联合优化<br/>最大化ELBO"] Loss2 --> Joint Joint --> Output["学习到的<br/>相互作用网络"] 相互作用的物理意义模型学习到的$K$种相互作用类型没有预先定义的物理含义，而是通过训练自动获得。通过对学习结果的分析，发现不同类型的相互作用对应不同的物理机制：强约束相互作用：对应于氢键、盐桥等强相互作用，限制残基相对运动弱耦合相互作用：对应于范德华力、疏水相互作用等弱相互作用，允许一定柔性动态介导相互作用：对应于在变构过程中变化的关键相互作用，如构象转换中的瞬时接触这种无监督学习方法避免了人为定义相互作用的局限性，能够发现传统方法难以识别的潜在相互作用模式。 Pin1系统：域间变构通信路径图2：Pin1在配体结合或突变时的蛋白质柔性和相互作用模式变化该图全面展示了Pin1在不同状态下的结构动力学和相互作用网络，是理解NRI模型如何从MD轨迹中学习变构信息的关键图示：图2a：蛋白质主链柔性变化（Backbone RMSD）具体内容：热图展示Pin1主链的均方根偏差（RMSD），颜色表示结构柔性颜色编码：蓝色（低RMSD，稳定）→红色（高RMSD，柔性）六种系统对比： apo-Pin1（无配体）：WW域（β1-β2）、催化环、α2螺旋和PPIase核心（β5/α4）显示高柔性（红色） FFpSPR-Pin1（正调控配体）：这些区域的柔性显著降低（变为蓝色），表明配体结合稳定了蛋白质构象 I28A突变：即使有FFpSPR结合，整体柔性增加，特别是WW域和催化环 pCdc25C-Pin1（负调控配体）：保持较高柔性，允许构象探索说明的问题：配体结合对柔性的影响：FFpSPR结合后，WW域和PPIase域的柔性被显著抑制正负调控差异：正调控配体使结构更刚性，负调控配体保持高柔性突变效应：I28A突变破坏了域间界面的稳定性逻辑链条：配体结合/突变 → 改变局部相互作用 → 影响结构柔性 → 反映在RMSD变化 → 指示变构效应存在图2b：残基间学习到的边缘分布图具体内容：点-线图，每个点代表一个残基，线代表NRI模型推断的显著相互作用表示方式：节点沿x轴排列，对应蛋白质序列位置边的颜色/粗细表示相互作用强度或类型说明的问题：相互作用网络拓扑：显示哪些残基对在动力学上耦合，即使它们空间距离可能较远 WW域的枢纽作用：WW域残基与其他区域有大量连接，表明其在动力学网络中的中心地位配体特异性模式：FFpSPR结合增强WW与PPIase核心间的连接，pCdc25C结合则产生不同的连接模式关键残基识别：I28、T29、C113等实验已知的重要位点在图中显示高连接度逻辑链条：NRI分析MD轨迹 → 推断残基间潜在相互作用 → 构建相互作用网络 → 识别网络中心和关键连接图2c：结构域/区块间边缘分布图具体内容：将相邻残基聚类为结构域/区块（如WW域、催化环、α1螺旋等），展示域间相互作用模式表示方式：矩阵热图或网络图，节点为结构域，边表示相互作用强度说明的问题：跨结构域通讯：显示哪些结构域在动力学上耦合，FFpSPR结合增强了WW与PPIase核心的连接变构通路可视化：清晰的域间连接模式，如WW→PPIase核心→催化环的路径调控机制差异：正调控增强域间连接，负调控减弱域间连接逻辑链条：残基水平相互作用 → 聚合到结构域水平 → 识别域间通讯模式 → 揭示变构调控的结构基础图2d：学习到的相互作用有向图具体内容：网络图表示，节点为结构域，边表示相互作用表示方式：节点大小：连接度（多少边连接到此节点）边粗细：相互作用强度箭头：影响方向（从发送方到接收方）说明的问题：信息流方向性：揭示变构信号的可能传递方向，如FFpSPR结合后信号从WW流向PPIase核心，再到催化环网络中心性分析：大节点是关键枢纽，如PPIase核心在多个系统中都是中心节点系统比较：不同配体/突变导致不同的网络拓扑，提供了变构机制的结构解释逻辑链条：NRI推断相互作用 → 构建有向网络 → 分析网络拓扑属性 → 推断信息流路径 → 解释变构机制综合逻辑链条整体分析框架：实验设计（不同配体/突变） MD模拟不同系统 NRI模型训练与推断相互作用图构建网络分析与通路识别机制解释与验证核心发现逻辑：变构信号传递路径的存在性证明：NRI成功推断出WW域到催化环的路径，这些路径在配体结合后增强，无配体时不存在正负调控机制对比：正调控（FFpSPR）增强域间连接，形成完整信号通路；负调控（pCdc25C）减弱域间连接，阻断信号传递突变效应解释：I28A突变破坏了WW与PPIase核心的连接，解释了其功能丧失方法优势验证：NRI能早期检测变构信号（50 ns内），比其他方法更敏感，能识别非线性、因果性相互作用 Pin1结构与功能 Pin1是一种包含两个结构域的肽酰脯氨酰顺反异构酶： WW域（残基1-39）：识别并结合磷酸化Ser/Thr-Pro基序，但无法催化异构化反应 PPIase域（残基50-163）：包含催化位点，执行肽酰脯氨酰键的顺反异构化 PPIase核心：α4-螺旋和β4-β7折叠片 α1-α3螺旋：形成催化位点的外壳催化环：半无序结构，参与底物结合和催化两个域通过连接肽（残基40-49）相连，形成独特的双域结构。WW域的结合能够变构调节PPIase域的活性，这种长程通信机制是Pin1功能调控的核心。配体结合的变构效应研究比较了五种状态的Pin1： apo-Pin1（PDB 3TDB）：无配体结合，WW域与PPIase域独立运动 FFpSPR-Pin1（PDB 3TDB）：正变构配体结合，WW域与PPIase域协调运动 I28A突变（PDB 3TDB）：域间界面突变，破坏WW-PPIase通信 pCdc25C-Pin1（PDB 1PIN）：负变构配体结合分离结构（PDB 1NMV）：WW域与PPIase域完全分离通过100 ns MD模拟（每20 ns采样一次，共50帧），NRI模型学习到了不同状态下的相互作用网络。关键发现： FFpSPR结合增强域间通信：学习到的边在WW域和其他结构域之间频繁出现，表明WW域是蛋白质运动的关键元素。具体表现为： WW域与PPIase核心之间的连接显著增强 WW域通过K97（α1-螺旋）和S105/C113（α2-3螺旋）与催化环建立新的通信路径域间界面（I28/T29）和催化位点附近（C113）的残基出现在变构路径上这些发现与实验研究一致，I28/T29和C113已被确定为影响Pin1活性的关键突变位点。图3：Pin1中介域间变构通信的路径通过计算学习到的网络中的最短路径，识别介导WW域到催化环的变构通信路径： (a) FFpSPR-Pin1的变构路径：三条路径从WW域出发，终结于催化环左侧路径：WW → Q131（PPIase核心）→ R69（催化环）中间路径：WW → P133（PPIase核心）→ S67（催化环）右侧路径：WW → K97（α1螺旋）→ S105/C113（α2-3螺旋）→ 催化环 (b) apo-Pin1：没有找到从WW域到催化环的路径，虽然WW域可以与α1-螺旋相互作用，但通信无法从α1-螺旋传递到催化环突变破坏域间通信 I28A突变的效应尤为显著：学习到的相互作用图显示，I28A突变急剧削弱了WW域与PPIase核心/α2-3螺旋之间的相互作用 WW域的涨落阻断了变构信号从WW向PPIase域的传播这表明I28在域间界面的关键作用，其突变导致蛋白质失去变构调控能力 pCdc25C结合的负变构效应： PPIase核心与WW域的相互作用减少 PPIase域内的边减少，反映域内接触减弱几乎没有边连接到催化环，表明PPIase域内的变构通信受阻分离结构（PDB 1NMV）的NRI分析：学习到的边主要集中在WW域与PPIase核心之间但与FFpSPR结合不同，WW域与α1-螺旋之间几乎无相互作用这表明空间接近但缺乏功能耦合时间依赖的信号传播通过分析不同时间窗口的相互作用演化，发现NRI模型能够在MD轨迹的早期阶段检测到变构信号： 50 ns（frames 1-500）：催化环中较大的边权重已被学习到 100 ns（frames 1-1000）：催化环的RMSD值增加3Å，反映连接到位点的边权重增强 200 ns（frames 1-2000）：传统的derivative centrality方法才能检测到完整的变构传播这表明NRI模型比传统方法提前数倍捕获变构信号，为理解变构机制提供了新的时间维度。 SOD1系统：突变诱导的构象变化图4：SOD1中G93A突变引起残基/域间相互作用变化该图揭示了与ALS相关的G93A突变如何通过变构机制影响SOD1的功能： (a) SOD1蛋白质的域划分：展示了G93A突变的位置（红色箭头）以及各个结构域 β桶（灰色）：8条反平行β折叠片，形成蛋白质核心二聚化环（DL，粉红色）二硫键环（DiL，绿色）锌结合环（ZL，橙色）静电环（EL，蓝色）：小的活性环 (b) WT SOD1和G93A SOD1在300 ns的初始结构： WT SOD1：EL稳定在金属位点附近（绿色箭头向上） G93A SOD1：EL远离金属位点（绿色箭头向下），表明构象变化 (c) WT（左）和G93A（右）在MD模拟中学习到的残基间边分布： WT：长活性环（DL、DiL、ZL）与小活性环（EL）紧密相互作用 G93A：长活性环内部连接几乎断裂，Zn（II）结合位点网络疏松 (d) 学习到的域间相互作用图： WT：活性环与β桶连接，导致EL闭合状态 G93A：活性环内连接断裂，EL开放 (e) 熵值归一化的边权重分布： WT：边权重集中在活性环内部 G93A：边权重分散，连接模式改变 (f) 从G93/A93开始的变构路径： WT（左）：G93 → DL → DiL → ZL → EL G93A（右）：A93 → β桶 → EL，不再通过长活性环 SOD1功能与ALS病理超氧化物歧化酶1（SOD1）是一种将超氧阴离子自由基转化为分子氧和过氧化氢的金属酶，在两步快速反应中交替还原和氧化活性位点铜。其整体结构由8条反平行β链加上形成活性位点的两个环组成。长活性环（残基49-83）可进一步分为：二聚化环（DL）：介导蛋白质二聚化二硫键环（DiL）：包含结构性二硫键锌结合环（ZL）：结合Zn（II）离子小活性环是静电环（EL），在金属位点附近发挥关键作用。 G93A突变与家族性肌萎缩侧索硬化症（ALS）相关：突变位点远离金属位点，属于典型的变构突变导致EL远离金属位点，降低Zn（II）亲和力影响ALS的病理过程 MD模拟与NRI分析对野生型（WT）和G93A SOD1进行500 ns MD模拟，分析结果：柔性变化： G93A SOD1的EL比WT更加柔性运动模式显示G93A突变诱导EL远离金属位点 WT SOD1的EL稳定在金属位点附近氢键网络： G93A突变使A93(O)-L38(N)距离增加，氢键相互作用减弱 β桶与活性环间的许多氢键被削弱 G93A SOD1结构比WT更加松散学习到的相互作用网络： WT SOD1：长活性环（DL、DiL、ZL）与小活性环（EL）紧密相互作用稳定Zn（II）结合环境长活性环和EL还连接到β桶中的残基，导致EL闭合状态变构路径从G93通过DL、DiL、ZL到EL G93A SOD1：长活性环内部的原始连接几乎断裂 Zn（II）结合位点网络疏松变构路径从A93直接通过β桶中的残基到EL，不再通过长活性环活性环内相互作用网络减弱，显著扩大Zn（II）结合口袋，降低Zn（II）亲和力这些发现完美解释了G93A突变的变构病理机制：通过破坏长活性环内的相互作用网络，导致Zn（II）结合环境不稳定，从而影响SOD1的催化功能和稳定性。 MEK1系统：激活相关的域通信 MEK1（MAPK/ERK激酶1）是RAS-RAF-MEK-ERK信号通路的关键组分，其活性受到多种机制的严格调控。研究了四种状态的MEK1： WT：野生型 A52V：非活性突变 E203K：活性突变（激活片段的螺旋-环转变） S218Sp/S222Sp：磷酸化激活（Ser218和Ser222磷酸化）通过MD模拟和NRI分析，揭示了激活相关的域间通信模式。结构域与激活机制 MEK1包含：小N叶：5条反平行β链（核心激酶域-1）和两个保守的αA/αC螺旋大C叶：3个核心激酶域、激活片段和富脯氨酸环激活片段的螺旋-环转变是MEK1激活的关键：非活性状态（WT、A52V）：激活片段为螺旋结构活性状态（E203K、S218Sp/S222Sp）：激活片段转变为环状结构学习到的相互作用网络 NRI模型揭示的域间通信模式：非活性MEK1（WT、A52V）：域间相互作用较少激活片段、富脯氨酸环与其他域的相互作用弱活性MEK1（E203K、S218Sp/S222Sp）： αA-螺旋、核心激酶域-1、激活片段和富脯氨酸环与其他域强烈相互作用这些域驱动磷酸化MEK1激活的慢速运动激活突变（E203K效应）：增强激活片段/富脯氨酸环与MEK1其他部分的相互作用从R201（近E203K）开始的变构路径显示，激活片段显著影响向富脯氨酸环传递信息通信通过αA-螺旋传播到αC-螺旋这些发现揭示了MEK1激活的变构机制：激活片段和富脯氨酸环形成相互作用模式，激活片段连接到αA-螺旋，可能影响其与激酶域其他部分的相互作用。方法优势与性能评估图7：基于Hessian和NRI的方法在捕获模拟中变构信号的性能对比该图对比了传统方法与NRI方法在检测变构信号方面的能力差异： (a, b) 基于Hessian的derivative node指标：在FFpSPR-和pCdc25C-Pin1系统中，使用轨迹不同片段计算δnode FFpSPR-Pin1：催化位点在200 ns（frame 2000）后才出现大的δnode值，表明完整的变构传播在200 ns后才被检测到 pCdc25C-Pin1：几乎没有信号传递到催化环，构象保持开放 (c, d) NRI方法学习到的域间边分布：显示域间相互作用和对应的平均构象（用RMSD值映射） FFpSPR-Pin1：50 ns（frames 1-500）内催化环中已学习到较大的边权重，开放构象在FFpSPR结合到WW域后约108 ns完成关闭转变 pCdc25C-Pin1：构象保持开放，几乎无信号传递到催化环早期信号检测 NRI模型的核心优势在于能够在MD轨迹的早期阶段检测到变构信号： 50 ns：NRI模型已在催化环中检测到较大的边权重 108 ns：开放构象完成关闭转变 200 ns：传统derivative centrality方法才检测到完整变构传播这表明NRI模型比传统方法提前约4倍时间捕获变构信号。自由能预测准确度图6：NRI方法计算自由能得分的性能评估该图验证了NRI方法在预测突变稳定性效应方面的准确性： (a) WT和23个Ala突变体的热力学数据总结，“N.D.”表示突变体太不稳定无法测量 (b) Ala突变对Pin1平衡稳定性的影响正值表示Ala突变相对于WT是去稳定的去稳定超过3 kcal/mol的突变显示为红色条，1-3 kcal/mol显示为蓝色条 (c, d) 基于NRI模型的计算自由能得分（ΔGZ）与实验自由能（ΔΔG）的对比 12Å相互作用阈值：$R^2 = 0.939$（95%置信区间：0.859 < $R^2$ < 0.974），$p = 3.361 \times 10^{-11}$ 15Å相互作用阈值：$R^2 = 0.931$（95%置信区间：0.842 < $R^2$ < 0.971），$p = 1.166 \times 10^{-10}$ (e) 基于约束网络分析（CNA）的计算自由能（ΔGCNA）与实验自由能的对比：$R^2 = 0.188$，$p = 0.390$ (f) MD模拟的总势能（ΔGTotal）与实验数据的对比：$R^2 = -0.093$，$p = 0.671$ 与传统方法的对比研究将NRI方法与三种传统方法进行了系统对比：方法原理局限性表现约束网络分析（CNA）基于Hessian的弹性网络模型假设设置，线性相关假设仅识别WW域的残基，遗漏催化环和α螺旋 Derivative centrality Hessian导数度量 200 ns后才检测到信号时间延迟显著动力学耦合指数（DCI）协方差矩阵替代Hessian 相关系数矩阵难以解读无法区分因果相关 NRI模型深度学习推断相互作用需要训练数据 50 ns检测信号，$R^2=0.939$ NRI模型的显著优势：早期检测：比传统方法提前数倍捕获变构信号因果推断：通过潜在变量建模相互作用，区分因果与非因果相关自由能预测：$R^2=0.939$ vs CNA的$R^2=0.188$，提升约5倍路径识别：能够识别多条变构路径，揭示冗余通信机制采样频率的影响研究系统评估了采样频率对学习结果的影响，使用10、15、20、25、30、40、50、60、75、90、100步进行测试：低频采样（≤50步）：产生相对较小的重建误差学习到的边较少且权重较低由于输入的结构信息较少，边的学习差异显著高频采样（>50步）：重建准确性显著下降采样间隔过大（如20步=250帧间隔）会错过许多关键的生物学功能构象权衡考虑：需要在采样频率和计算效率之间权衡步长间隔约20 ns可产生更合理的结果基于小的重建误差和充分采样选择学习结果模型消融实验为测试图神经网络在NRI中的作用，进行了消融实验，将提出模型与无潜在边变量的变分自编码器（VAE）基线进行对比：将轨迹分割为训练/验证/测试集 Pin1、MEK1和SOD1的MSE结果显示，边上的潜在变量改善了模型性能提出的架构为MD轨迹的边（残基相互作用）建模提供了更好的框架在密集相互作用系统中（如WT-SOD1），NRI模型的优势更加显著 Q&A Q1：NRI模型与传统MD分析方法（如RMSD、RMSF、PCA）有什么本质区别？为什么深度学习方法能捕获传统方法难以识别的信息？ NRI模型与传统MD分析方法的根本区别在于信息提取方式和因果推断能力：分析方法提取信息局限性适用场景 RMSD/RMSF 整体/局部结构变化无法区分长程通信，忽略因果判断平衡、识别柔性区域 PCA/EFA 主要运动模式线性组合，难以捕获非线性相互作用构象态聚类互相关分析残基间相关性无法区分因果vs非因果相关初步识别关联 NRI模型因果相互作用网络需要训练数据识别变构路径、预测自由能深度学习的独特优势：非线性建模能力：NRI通过GNN的message passing机制，能够捕获残基间复杂的非线性相互作用，而传统方法通常基于线性假设或弹性网络模型。因果推断：NRI通过潜在变量$z$建模相互作用，并通过重建任务验证其有效性。这确保学习到的是对系统演化有因果贡献的相互作用，而非仅仅是统计相关。高维特征抽象：NRI的encoder将高维轨迹（$3N$维）映射到低维潜在空间（$K$种相互作用类型），自动提取对系统演化最关键的特征。动态网络视角：将蛋白质变构建模为动态演化的相互作用网络，而非静态结构或单一势能面，更符合生物系统的本质。形象类比：传统方法：像是拍摄交通视频后统计每辆车的速度和位置，但无法识别“交通瓶颈” NRI模型：像是分析车与车之间的相互作用（跟车、变道、超车），识别出“一旦堵塞就会导致全城瘫痪”的关键路口（变构热点） Q2：NRI模型学习到的K种相互作用类型是否有明确的物理意义？如何解释不同类型的相互作用？ NRI模型学习到的$K$种相互作用类型没有预先定义的物理含义，但通过训练自动获得了明确的物理意义。这是一种无监督学习的优势：避免了人为定义相互作用的偏差和局限性。相互作用类型的物理意义通过对三个系统（Pin1、SOD1、MEK1）学习结果的分析，可以归纳出以下几种典型的相互作用类型：相互作用类型物理意义特征出现位置强约束型氢键、盐桥、π-π堆积边权重大，在所有状态下稳定二级结构内部、结构域核心弱耦合型范德华力、疏水相互作用边权重小，波动较大结构域界面、loop区动态介导型变构过程中瞬时接触仅在特定状态出现变构路径上稳定抑制型空间位阻、排斥作用负边权重，减少运动构象转换的屏障协同增强型别构效应增强边权重随时间增加配体结合后的域间通信在Pin1系统中的具体体现在FFpSPR-Pin1的NRI分析中，观察到的相互作用类型模式：类型1-3：在WW域和PPIase核心之间的高权重边物理意义：域间界面的氢键网络和疏水核心功能：稳定双域结构，介导长程通信类型4-6：在α1/α2-3螺旋与催化环之间的中等权重边物理意义：变构通信的关键桥梁功能：传递信号从WW域到催化位点类型7-10：在PPIase域内部的低权重边物理意义：柔性调节和构象涨落功能：允许必要的构象变化在SOD1系统中的具体体现在WT vs G93A SOD1对比中，相互作用类型的显著差异： WT SOD1：类型1-4主导：长活性环（DL、DiL、ZL）内部强相互作用物理意义：稳定Zn（II）结合环境功能：维持EL闭合状态 G93A SOD1：类型5-8出现：β桶与EL之间的直接相互作用类型1-4显著减弱：长活性环内部连接断裂物理意义：变构突变导致相互作用网络重排功能：导致EL开放，Zn（II）亲和力降低验证相互作用类型的有效性通过以下方式验证学习到的相互作用类型的物理意义：与已知实验数据对比：学习到的关键残基（如Pin1的I28/T29/C113）与实验验证的变构热点一致自由能预测准确度：基于学习到的相互作用网络计算的自由能变化与实验数据高度相关（$R^2=0.939$）时间一致性检验：在重复的MD模拟中，学习到的相互作用拓扑高度一致，特别是关键的拓扑元素（如MEK1的激活片段和富脯氨酸环）消融实验：移除边潜在变量后的VAE基线模型性能下降，证明边上的潜在变量捕获了真实的物理相互作用未来改进方向虽然NRI模型学习到的相互作用类型具有明确的物理意义，但可以通过以下方式进一步增强可解释性：有监督训练：使用已知的相互作用类型（如氢键、盐桥）作为标签，使模型直接学习这些类型后验分析：对每个相互作用类型的残基对进行结构分析，归纳共同的几何和物理化学特征注意力机制：在GNN中引入注意力权重，提供更细粒度的相互作用强度解释 Q3：NRI模型对采样频率和轨迹长度有什么要求？如何确定合适的采样参数？ NRI模型对采样频率和轨迹长度的要求需要仔细权衡，这涉及MD模拟的计算成本和模型学习效果的平衡。采样频率的影响研究系统测试了10、15、20、25、30、40、50、60、75、90、100步的采样间隔，发现了以下规律：低频采样（≤50步）：优势：重建误差（MSE）和方差相似度（VSD）较小计算效率高劣势：学习到的边较少且权重较低由于输入结构信息较少，边的学习差异显著对于构象变化显著的系统（如pCdc25C-Pin1），学习结果不稳定高频采样（>50步）：优势：输入信息更丰富学习结果更稳定劣势：重建准确性显著下降采样间隔过大可能错过关键构象计算成本高临界阈值：采样间隔约20 ns是一个合理的上限超过20 ns可能太长，无法恢复变构过程中的足够信息例如，选择20步会导致250帧的间隔，错过许多关键的生物学功能构象推荐的采样策略基于研究结果，推荐以下采样策略：系统类型推荐采样间隔轨迹长度采样帧数理由快速变构系统（如Pin1） 10-20 ns 100-200 ns 10-20帧捕获快速构象转变慢速变构系统（如SOD1） 20-40 ns 500 ns 15-25帧平衡采样密度和计算成本突变效应研究 20 ns 200-500 ns 10-25帧捕获突变前后差异轨迹长度的影响研究对不同时间窗口的边分布进行了分析：滑动窗口分析（frames 1-1000, 1000-2000, …, 4000-5000）：生物分子的动力学随时间显著变化不同时间段的边分布差异较大累积窗口分析（frames 1-500, 1-1000, …, 1-5000）：边分布相对稳定反映整个动态过程的整体特征，而非每个片段的特征推荐策略：使用累积窗口（frames 1-N）进行分析确保轨迹长度足够捕获至少一次完整的构象转变对于Pin1，100-200 ns足够捕获open-to-closed转变对于SOD1，500 ns足够捕获突变诱导的构象变化模型训练的稳定性研究进行了三次重复MD模拟，验证了NRI模型的稳定性： Pin1系统：重复轨迹的边分布相似但有差异基础拓扑（WW→PPIase核心）稳定 SOD1系统：重复轨迹的边显示高度一致性表明NRI模型在WT-SOD1情况下捕获边更准确 MEK1系统：边的差异略大但重要的拓扑元素（激活片段和富脯氨酸环）学习一致实际应用建议基于研究结果，实际应用NRI模型的建议：初步探索：使用较短轨迹（100-200 ns）和较高采样频率（10-20 ns）快速评估系统的变构行为精细分析：使用较长轨迹（500 ns）和中等采样频率（20-40 ns）平衡计算成本和学习效果验证策略：检查VSD值，确保重建误差可接受（VSD < 0.2）进行重复模拟，验证学习结果的稳定性对比不同采样间隔的结果，选择最优参数计算资源有限时：优先保证采样频率而非轨迹长度过长的低频采样轨迹可能不如适中的高频采样轨迹关键结论与批判性总结核心贡献深度学习赋能MD分析：首次将神经关系推断（NRI）模型应用于生物分子MD数据分析，通过图神经网络同时推断残基间的潜在相互作用，将蛋白质变构过程建模为动态演化的相互作用网络早期信号捕获：NRI模型能够在MD轨迹的早期阶段（50-100 ns）检测到变构信号，比传统基于Hessian的方法（200 ns以后）提前数倍，为理解变构机制提供了新的时间维度自由能准确预测：基于学习到的相互作用网络计算突变后的相对自由能变化，与实验数据高度一致（$R^2=0.939$，$p=3.361 \times 10^{-11}$），显著优于传统约束网络分析（CNA）方法（$R^2=0.188$，$p=0.390$）多系统验证：在Pin1（域间变构）、SOD1（突变病理）、MEK1（激活机制）三个不同的变构系统中成功识别长程通信路径，证明了方法的普适性物理可解释性：学习到的相互作用类型具有明确的物理意义（强约束、弱耦合、动态介导等），能够识别实验验证的关键残基（如Pin1的I28/T29/C113）局限性采样频率敏感性：NRI模型对采样频率较为敏感，低频采样（≤50步）虽然计算效率高但可能遗漏关键构象，高频采样（>50步）计算成本高且重建误差大。需要根据具体系统在采样密度和计算效率之间权衡轨迹长度要求：虽然NRI能在早期阶段检测到变构信号，但仍需要足够长的轨迹（100-500 ns）来捕获完整的构象转变和达到统计收敛。对于慢速变构系统（毫秒级），常规MD仍无法覆盖完整过程因果推断的隐含假设：NRI通过重建任务验证相互作用的有效性，但重建误差小不一定等同于因果关系的正确性。可能存在一些在重建任务中不重要但在生物学功能上关键的相互作用被遗漏黑箱模型的解释性：虽然学习到的相互作用类型具有物理意义，但GNN的decision-making过程仍是黑箱，难以完全解释为何特定残基对被归类为某种相互作用类型超参数选择：模型包含多个超参数（相互作用类型数$K$、GNN层数、隐藏维度等），文中未详细讨论这些参数的选择原则和对结果的影响未来研究方向扩展到更大尺度系统：研究NRI模型在多亚基蛋白复合物、蛋白质-核酸复合物、超大分子组装体（如核糖体、蛋白酶体）中的表现，评估其在更复杂系统中的泛化能力整合多尺度建模：结合增强采样技术（如加速MD、Metadynamics）或马尔可夫态模型（MSM），将NRI的应用范围扩展到毫秒-秒级的慢速变构过程有监督相互作用分类：使用已知的相互作用类型（氢键、盐桥、π-π堆积等）作为标签，使模型直接学习这些类型，进一步增强可解释性实时变构监测：开发在线学习版本的NRI，能够在MD模拟过程中实时更新相互作用网络，实现变构信号的实时监测和预警结合实验数据：整合NMR、HDX-MS、FRET等实验数据作为约束或验证，提高学习到的相互作用网络的准确性和生物学相关性方法比较与基准测试：在更多蛋白质家族和变构类型中系统比较NRI与其他深度学习方法（如VAE、GAN、Transformer），建立标准化的评估基准药物设计应用：将NRI识别的变构热点和通信路径用于变构药物设计，预测和优化变构调节剂的结合位点代码与工具开发：虽然论文提供了GitHub代码，但需要进一步开发用户友好的软件包和可视化工具，降低方法使用门槛，使更多研究者能够应用NRI解决实际问题小编锐评：这篇文章的核心思路很清晰：用NRI把MD轨迹变成相互作用网络，然后从中挖掘变构路径和自由能变化最吸引人的是能在50-100 ns检测到变构信号，比传统方法快4倍，这对MD模拟来说意义重大但文章对模型超参数选择、不同深度学习架构的系统比较讨论较少，是未来研究可以补充的地方 $R^2=0.939$的自由能预测确实很惊艳，但只在Pin1的23个Ala突变上验证，还需要在更多系统上测试代码开源了，但不知道易用性如何，希望有更友好的界面让非计算机背景的研究者也能用
Molecular Dynamics · 2026-01-25
变构激活的动态基础：恶性疟原虫蛋白激酶G的长程通信机制
 变构激活的动态基础：恶性疟原虫蛋白激酶G的长程通信机制本文信息标题: 变构激活的动态基础：恶性疟原虫蛋白激酶G的长程通信机制作者: Jinfeng Huang, Jung Ah Byun, Bryan VanSchouwen, Philipp Henning, Friedrich W. Herberg, Choel Kim, Giuseppe Melacini 发表时间: 2021年6月10日单位: McMaster University（加拿大麦克马斯特大学）, University of Kiel（德国基尔大学）, Baylor College of Medicine（美国贝勒医学院）, Rice University（美国莱斯大学）引用格式: Huang, J., Byun, J. A., VanSchouwen, B., Henning, P., Herberg, F. W., Kim, C., & Melacini, G. (2021). Dynamical Basis of Allosteric Activation for the Plasmodium falciparum Protein Kinase G. The Journal of Physical Chemistry B, 125(23), 6532-6542. https://doi.org/10.1021/acs.jpcb.1c03622 摘要恶性疟原虫的cGMP依赖性蛋白激酶（PfPKG）对于疟原虫生命周期的进程是必需的，因此是一个有前景的抗疟药物靶点。PfPKG包含四个cGMP结合结构域（CBD-A至CBD-D）。CBD-D在PfPKG调控中发挥关键作用，它是催化结构域抑制和cGMP依赖性激活的主要决定因素。因此，理解CBD-D如何被cGMP变构调节至关重要。虽然CBD-D的apo与holo构象变化已有报道，但目前缺乏关于激活途径中间态的信息。在本研究中，我们采用分子动力学模拟来建模PfPKG CBD-D结构域cGMP依赖性激活热力学循环中的四个关键状态。模拟结果与NMR数据进行比较，揭示了PfPKG CBD-D激活途径会采样一种紧凑中间态，其中N端和C端螺旋靠近中央β桶。此外，通过比较cGMP结合的活性态和非活性态，识别了区分这两种状态的关键结合相互作用。识别cGMP结合非活性态特有的结构和动力学特征，为设计PfPKG选择性变构抑制剂作为疟疾的可行治疗方案提供了有希望的基础。核心结论四态热力学循环：首次完整映射了PfPKG CBD-D的变构激活路径，包括难以捕捉的apo/active和holo/inactive中间态区域特异性响应：PBC区域的动力学抑制需要cGMP结合和变构构象变化的协同作用，而αB-αC螺旋主要由变构效应调控变构抑制剂设计基础：holo/inactive中间态的结构特征，特别是R484-A485与cGMP相互作用的变化，为设计选择性变构抑制剂提供了明确靶点物种选择性机制：PfPKG的R484与人类PKG的K308在αC螺旋相互作用上的差异，可实现宿主-寄生虫选择性背景关键术语解释在深入讨论之前，先介绍本文涉及的关键缩写： PfPKG：Plasmodium falciparum cGMP-dependent protein kinase G（恶性疟原虫cGMP依赖性蛋白激酶G） cGMP：cyclic guanosine monophosphate（环磷酸鸟苷），细胞内第二信使分子 CBD：cGMP-binding domain（cGMP结合结构域），负责识别和结合cGMP PBC：Phosphate-Binding Cassette（磷酸结合盒），CBD中结合cGMP磷酸基团的关键区域 BBR：Base-Binding Region（碱基结合区），CBD中结合cGMP鸟嘌呤碱基的区域 N3A：N-terminal three-helix assembly（N端三螺旋组装体），包含αX:N、α310和αA螺旋的复合结构 apo：配体未结合状态（如无cGMP结合的蛋白状态） holo：配体结合状态（如cGMP结合的蛋白状态） β-core：中央β桶，CBD结构域的核心支架，由8个β折叠片组成 cation-π相互作用：阳离子-π相互作用，带正电荷的离子（如铵根离子）与芳香环的π电子云之间的静电相互作用，在蛋白质-配体识别中很重要 His τ态中性：组氨酸在pH=7时的质子化状态，质子位于Nε2（τ氮）上，整体不带电（记为HIE），是生理条件下最常见的组氨酸状态，适用于大多数蛋白质MD模拟疟疾与PfPKG的重要性疟疾是由恶性疟原虫（Plasmodium falciparum）引起的致命寄生虫病，每年导致全球数十万人死亡。疟原虫的生命周期复杂，包括在蚊虫中的有性生殖阶段和在人体内的无性增殖阶段，其中从肝细胞释放出的裂殖子侵入红细胞是引发疟疾症状的关键步骤。 PfPKG是一个cGMP依赖性丝氨酸/苏氨酸激酶，在疟原虫的生命周期调控中扮演中央开关的角色。研究表明，PfPKG在疟原虫的多个关键生命周期阶段都发挥着不可替代的作用，包括裂殖子从红细胞释放（egress）、裂殖子重新侵入红细胞（invasion）以及配子体激活（sexual stage development）。抑制PfPKG的活性可以阻断这些关键过程，从而阻止疟原虫的生命周期进程，因此PfPKG被认为是极具前景的抗疟药物靶点。特别值得注意的是，PfPKG与人类PKG在结构上存在差异，这为实现宿主-寄生虫选择性抑制提供了可能性，即可以设计只杀灭疟原虫而不伤害人体正常细胞的药物。 cGMP结构域与变构激活机制 PfPKG包含四个cGMP结合结构域（CBD-A、CBD-B、CBD-C和CBD-D），位于N端调控区，其中CBD-D具有最高的cGMP结合亲和力（Kd = 51 ± 7 nM），是变构调控的核心决定因素。此外，PfPKG还包含一个催化结构域，位于C端，负责ATP（Adenosine Triphosphate，三磷酸腺苷，细胞能量货币和磷酸供体）结合和磷酸转移反应，在无cGMP状态下被N端结构域抑制，cGMP结合后解除抑制。在无cGMP状态下，CBD结构域与催化结构域通过αB-螺旋和连接区相互作用，抑制催化活性。当cGMP结合到CBD-A和CBD-B时，引发变构激活：CBD-A结合cGMP解除对催化结构域的抑制，而CBD-B结合cGMP进一步激活催化结构域。然而，这一过程的原子级动态机制和长程通信路径尚未明确，尤其是连接apo/inactive到holo/active转变的中间态（如apo/active和holo/inactive）仍难以通过实验手段表征。变构激活的科学问题经典变构理论认为，配体结合通常稳定化蛋白局部结构，从而引发下游效应。但对于PfPKG，存在多个尚未解决的关键问题：CBD-A和CBD-B的cGMP结合是否都导致局部稳定化，还是存在区域特异性差异？局部变化如何跨越约60Å的距离传播至催化结构域，具体的信号传播路径是什么？催化结构域的哪些区域对变构信号最敏感，这些区域的动态变化如何与激酶活性相关？这些问题需要结合实验动态测量（如NMR化学位移分析）和原子级模拟（如微秒级MD模拟）来回答，特别是需要表征难以捕捉的中间态（如apo/active和holo/inactive）。关键科学问题本研究重点关注三个关键科学问题。四态变构循环的动态特征问题涉及PfPKG CBD-D的激活途径是否遵循离散的四态模型（apo/inactive、apo/active、holo/inactive、holo/active），以及不同状态间的转变路径和能量景观如何分布。区域特异性的变构响应问题关注PBC和αB-αC螺旋对cGMP结合和变构效应的敏感性是否存在显著差异，以及这种差异如何影响变构信号传播。变构抑制剂的设计基础问题则探索holo/inactive中间态具有哪些独特的结构和动力学特征，以及如何利用这些特征设计可结合但不激活激酶的选择性变构抑制剂，同时实现对PfPKG和人类PKG的区分。创新点方法学创新：首次将NMR实验与MD模拟结合研究PfPKG完整四态变构循环，实验-计算互补验证动态变化中间态表征：首次在原子分辨率下表征了难以捕捉的apo/active和holo/inactive中间态变构抑制剂设计基础：识别了holo/inactive中间态的独特结构特征，为设计可结合但不激活的选择性抑制剂提供了明确靶点区域特异性机制：揭示了PBC和αB-αC螺旋对cGMP结合和变构效应的不同敏感性，深化了对变构通信机制的理解图S1：四态变构循环的初始结构模型四态初始结构的建模本研究仅两态有实验解析的晶体结构，另外两态通过计算建模获得：实验解析的晶体结构 apo/inactive状态：PDB 4OFF（apo CBD-D晶体结构） holo/active状态：PDB 4OFG（cGMP-bound CBD-D晶体结构）计算建模的中间态状态建模方法结构来源关键操作 apo/active 从holo/active移除cGMP 4OFG 移除cGMP，保留活性构象（N3Aout/BCin） holo/inactive cGMP对齐到inactive结构 4OFF + 4OFG 通过β-core区域对齐，将cGMP从4OFG对齐到4OFF apo/inactive (补充) 添加缺失残基 4OFF + 5DYK 从全长结构(PDB 5DYK)补充N端2个残基和C端残基517-542 关键建模细节 apo/active状态：直接从holo/active晶体结构（4OFG）中移除cGMP，保持活性构象（N3Aout/BCin拓扑） holo/inactive状态：将holo/active（4OFG）和apo/inactive（4OFF）结构在保守的β-core区域对齐，然后将4OFG中的cGMP分子转移到4OFF结构中，创建一个配体结合但不激活的模型 apo/inactive补充：4OFF结构缺失N端前2个残基和C端517-542残基，从全长apo/inactive结构（PDB 5DYK）移植这些缺失区域，并通过β-core对齐确保结构连续性这种建模策略使得MD模拟能够探索难以通过实验表征的中间态（apo/active和holo/inactive），从而完整映射四态变构热力学循环。研究方法：NMR与MD模拟的结合本研究采用实验-计算双管齐下的策略：核磁共振（NMR）实验测量野生型和突变型PfPKG CBD-D在cGMP结合状态下的化学位移通过化学位移导出的序参量（$S^2$，Order Parameter）评估蛋白质骨架动力学，$S^2$值范围0-1，越接近1表示运动越受限比较不同变构状态下的NMR数据，识别关键构象变化突变实验验证MD模拟预测的关键相互作用图S2：MD模拟与NMR实验的验证对比了三种力场（FF99SBnmr、FF14SB、FF99SBildn）预测的N-H序参量（$S^2$）与NMR实验数据黑色点为NMR实验值，绿色/红色/蓝色条为不同力场的MD预测值垂直箭头标注实验观察到的局部极小值结论：FF99SBnmr力场与实验数据最为一致，因此作为后续分析的主力场分子动力学（MD）模拟对四态变构循环中的每个状态进行3×1 μs重复模拟（总计12 μs）分析均方根偏差（RMSD，Root Mean Square Deviation），衡量结构与参考构象的偏离程度分析均方根涨落（RMSF，Root Mean Square Fluctuation），衡量原子运动的柔性使用CHESPA（Chemical Shift Projection Analysis，化学位移投影分析）比较突变效应通过相似性测量（SM，Similarity Measure）图谱映射构象转变路径 MD模拟细节使用Amber 16与GPU版pmemd.cuda在SHARCNET平台运行 cGMP参数通过HF/6-31G*量子化学计算获得电荷，经RESP（Restrained Electrostatic Potential，限制静电势）拟合得到部分电荷，并采用GAFF（General Amber Force Field，通用AMBER力场）补全缺失参数蛋白使用FF99SBnmr（专门为NMR数据优化的AMBER力场）为主力场，FF99SBildn（改进的侧链二面角参数）与FF14SB（AMBER 2014力场）用于holo/active对照体系溶剂化于TIP3P水盒子，边界距溶质至少12 Å；加入NaCl至100 mM模拟生理盐浓度 pH设为7，His为τ态中性（质子位于Nε2，记为HIE）；N/C端与Asp/Glu/Arg/Lys为标准电离态四态构象各进行3×1 μs轨迹，另对holo/active用两种力场各补充3 μs，总计18 μs 能量最小化后分段升温与平衡：NVT 0–100 K（20 ps），NPT 100–306 K（80 ps），逐步降低主链约束生产期在306 K、1 atm的NPT条件下运行，非键截断12 Å，长程静电相互作用用PME（Particle Mesh Ewald，粒子网格Ewald方法）轨迹每10 ps存储一次，分析使用CPPTRAJ（Amber工具包中的轨迹分析程序）结果与讨论 1. CBD-D结构域的动态分析图2：PfPKG CBD-D四态的全蛋白主链RMSD随时间变化 (A-D) 四态的RMSD时间轨迹：(A) Apo/Inactive，(B) Apo/Active，(C) Holo/Inactive，(D) Holo/Active 计算方法：将整个蛋白的主链（N、Cα、C原子）对齐到各自状态的初始模型，计算RMSD 横轴为模拟时间（ns），纵轴为RMSD（Å）每个状态有3条1 μs独立轨迹，用不同灰度表示（黑色、深灰、浅灰）关键发现：所有12条轨迹（四态×3次重复）在1 μs内保持稳定，没有持续上升或大的构象漂移，表明模拟已达到平衡，可用于后续分析 RMSF：残基级别的柔性变化均方根涨落（RMSF）分析揭示了四态变构循环中的区域特异性动态响应。通过overlay整个CBD-D的Cα原子到初始模型，计算每个残基的RMSF值，发现：图3：PfPKG CBD-D残基特异性结构涨落（RMSF） (A) 全域RMSF vs 残基编号，四态用不同颜色表示：红色（apo/inactive）、蓝色（apo/active）、橙色（holo/inactive）、绿色（holo/active）。灰色高亮显示四态间最显著差异的区域，y轴使用log10刻度 (B-E) 不同状态对间的RMSF差异图：B和C量化变构构象变化的效应，D和E量化cGMP结合的效应关键发现：PBC和αB-αC螺旋对变构信号和cGMP结合的敏感性截然不同区域特异性RMSD分布为进一步量化不同结构元件的动态变化，研究分别计算了N3A区域、PBC区域和αB-αC螺旋的RMSD分布（通过overlay各自的β-core到初始结构，确保仅测量局部构象变化）。图4：N3A、PBC与αB-αC区域的特异性动态响应 (A-C) 分别展示N3A、PBC、αB-αC区域的RMSD箱线图，通过overlay β-core到初始模型计算。横轴为四态，纵轴为RMSD（Å） (D) 全域RMSD分布（overlay整个CBD-D主链到初始结构）箱线图说明：中线为中位数，箱体为25%-75%分位数，须为1.5×IQR范围，小方块为均值，两个叉号为1%和99%分位数区域四态RMSD特征调控机制物理意义 N3A (图4A) 四态间分布相似由整体构象决定，而非cGMP结合 N3A的in/out取向在所有状态下都能动态采样，与β-core的相对位置稳定 PBC (图4B) holo/active显著低于其他三态 cGMP结合和变构激活的协同作用 PBC稳定化需要双重因素，验证了RMSF结果 αB-αC螺旋 (图4C) active状态低于inactive状态主要由变构效应决定 αB-αC螺旋的动态性主要受构象状态调控，cGMP结合影响较小全域 (图4D) 反映αB-αC的大幅变化变构贡献占主导因αB-αC构象变化幅度最大，全域RMSD主要反映其变化 2. 变构转变路径：从inactive到active SM图谱的计算方法相似性测量（SM，Similarity Measure）是一种基于RMSD的二维散点图，用于直观评估构象在active和inactive状态之间的相对位置。对MD轨迹中的每一帧构象，分别计算： \[X = \mathrm{RMSD}_{\mathrm{N3A}}^{\mathrm{active}} - \mathrm{RMSD}_{\mathrm{N3A}}^{\mathrm{inactive}} \\ Y = \mathrm{RMSD}_{\alpha\mathrm{B}\text{-}\alpha\mathrm{C}}^{\mathrm{active}} - \mathrm{RMSD}_{\alpha\mathrm{B}\text{-}\alpha\mathrm{C}}^{\mathrm{inactive}}\] 符号区域相对于谁的RMSD 参考结构 $\mathrm{RMSD}_{\mathrm{N3A}}^{\mathrm{active}}$ N3A区域 active结构 holo/active晶体（PDB 4OFG） $\mathrm{RMSD}_{\mathrm{N3A}}^{\mathrm{inactive}}$ N3A区域 inactive结构 apo/inactive晶体（PDB 4OFF） $\mathrm{RMSD}_{\alpha\mathrm{B}\text{-}\alpha\mathrm{C}}^{\mathrm{active}}$ αB-αC螺旋 active结构 holo/active晶体（PDB 4OFG） $\mathrm{RMSD}_{\alpha\mathrm{B}\text{-}\alpha\mathrm{C}}^{\mathrm{inactive}}$ αB-αC螺旋 inactive结构 apo/inactive晶体（PDB 4OFF）计算步骤：对MD轨迹的每一帧，分别计算N3A和αB-αC区域相对于active和inactive参考结构的RMSD 计算差值得到 $(X, Y)$ 坐标在二维平面上绘制每帧的坐标点图5：PfPKG CBD-D的活性-非活性转变路径映射 (A, B) N3A与αB-αC的RMSD相似性测量（SM）图谱，展示apo/inactive（红色）和apo/active（蓝色）模拟轨迹。每个象限代表N3A和αB-αC结构元件的不同in/out组合姿态。A和B面板仅在数据集的前后显示顺序上不同 (C, D) 与A、B相同，但展示holo/inactive（橙色）和holo/active（绿色）模拟轨迹 (E) 总结PfPKG CBD-D沿变构热力学循环的主要动态变化的示意图。实线表示inactive（红色）和active（绿色）状态的初始拓扑结构，虚线和黑色箭头表示转变过程中的主要拓扑变化这种作差的方法勉强可借鉴吧，甚至可以作为CV？这种模拟也算是类似于，用增强采样采到了一些关键态，再跑standard MD得到kinetics 象限映射与物理意义象限坐标构象组合代表的状态拓扑特征右上 (+, +) N3Aout/BCin Holo/active参考态 N3A向外，αB-αC向内（活性）左下 (-, -) N3Ain/BCout Apo/inactive参考态 N3A向内，αB-αC向外右下 (+, -) N3Ain/BCin 紧凑中间态两者都向内，过渡态的必经之路（最多采样）左上 (-, +) N3Aout/BCout 松散中间态两者都向外（较少采样） Figure 5的SM图谱揭示了PfPKG CBD-D变构激活的能量景观。四个象限代表四个不同的构象 basin，每个数据点代表MD轨迹中的一帧构象。象限偏好性反映能垒： apo/inactive轨迹（红色）：主要分布在左下象限（N3Ain/BCout），与初始构象一致，表示inactive状态是稳定的能量极小值 holo/active轨迹（绿色）：主要分布在右上象限（N3Aout/BCin）和右下象限，表明active状态虽以N3Aout/BCin为主，但会大量采样紧凑中间态紧凑中间态的关键作用：右下象限（N3Ain/BCin）的数据点密度最高，所有四态的轨迹都显示出对这个象限的偏好采样这个紧凑中间态是inactive→active转变的必经之路，在能量景观中代表一个能量较低的区域物理上，N3Ain/BCin构象具有最小的空间位阻，是结构重排的最优路径松散中间态的稀有性：左上象限（N3Aout/BCout）的采样最少，表明松散构象在能量上不利这可能是因为N3Aout/BCout构象导致空间位阻增大，或者破坏了关键的分子内相互作用与PBC视角的一致性（Figure S3）：当用PBC替换N3A进行SM分析时（Figure S3），观察到相似的象限偏好性：所有激活路径都偏好紧凑的PBCin/BCin中间态（注意：PBC的in对应active构象），而非松散的PBCout/BCout路径。这进一步验证了紧凑中间态的普适性。图S3：PBC视角的活化-非活化转变路径 (A-B) Apo状态的PBC vs αB-αC SM图谱，比较PBC与αB-αC区域在active与inactive结构间的差异 (C-D) Holo状态的SM图谱，展示相同区域的构象变化关键发现：与Figure 5类似，所有激活路径都偏好紧凑的PBCin/BCin中间态，而非松散的PBCout/BCout路径重要结论基于Figure 5和S3的SM图谱分析，我们得出以下关键结论：紧凑中间态是变构转变的瓶颈：Figure 5的SM图谱显示所有四态轨迹都对右下象限（N3Ain/BCin紧凑中间态）有偏好采样，数据点密度最高。文献基于此推论认为这是inactive→active转变的”obligatory”（必经）中间态，物理上对应最小的空间位阻。需要注意的是，SM图谱本身不能直接观察完整的转变路径，这一推论仍需单分子实验或毫秒级增强采样进一步验证。 apo/active中间态的混合特征：结合了holo/active和apo/inactive的元素——PBC动力学类似apo/inactive（较不稳定，需要cGMP结合来稳定），而αB-αC螺旋构象类似holo/active（较稳定，主要由变构状态调控）。这解释了为什么apo/active状态的SM分布跨越多个象限。 holo/inactive中间态更接近inactive：无论在PBC还是αB-αC水平，holo/inactive都更像apo/inactive而非holo/active。这表明单靠cGMP结合不足以驱动active构象，必须同时满足变构构象变化才能实现激活，验证了PBC的双重依赖机制。 N3A的动态采样特性：N3A在所有四个状态下都能动态采样in和out取向（Figure 5E显示N3A的双向箭头），这与其在结构上的相对独立性有关。相比之下，αB-αC螺旋的in/out转变更受构象状态约束（Figure 4C显示active状态αB-αC更稳定）。 3. C端螺旋相互作用：激酶激活的关键接触与人类PKG和HCN通道的比较图S5：PfPKG与人类PKG的αC螺旋相互作用对比 (A) Holo/Active的PfPKG CBD-D（N3Aout/BCin）与人类PKG Iβ CBD-B的叠合视图。PfPKG用绿色丝带表示，人类PKG Iβ用青色丝带表示，cGMP与关键残基以棒状显示。两者在β-core上对齐，便于比较lid区域与αC螺旋的接触 (B) Holo/Inactive的PfPKG CBD-D（N3Ain/BCout）与人类PKG Iβ CBD-B的叠合视图。PfPKG以橙色系表示，人类PKG Iβ以浅色半透明丝带表示，cGMP与关键残基以棒状显示，用于对比非活化构象下的lid位置与cGMP周围相互作用关键差异：PfPKG的R484可与C端αC螺旋Q532/D533形成capping triad，而人类PKG Iβ对应的K308不形成类似稳定接触，为选择性变构抑制提供了结构依据两个面板均以β-core为对齐基准，强调lid与αC螺旋相互作用的物种差异 PfPKG的变构机制与哺乳动物PKG存在显著差异。人类PKG Iβ的CBD-B中，αB-螺旋在cGMP结合后动力学降低（保护作用），而PfPKG的CBD-B显示动力学增强（去保护作用）。这种差异使得CBD-B成为PfPKG选择性抑制的潜在靶点。与HCN（超极化激活环核苷酸门控）通道相比，PfPKG的变构转变路径更为单一，所有激活路径都经过“紧凑”N3Ain/BCin中间态，而HCN遵循多分支的路径。这表明不同环核苷酸结合结构域的变构调控机制存在显著多样性。关键相互作用通过比较holo/active和holo/inactive状态的N3Aout/BCin和N3Ain/BCout构象，可以识别激酶激活所需的关键相互作用。图6：C端螺旋与PBC的相互作用分析 (A, E) PfPKG CBD-D C端αC螺旋与PBC、Y480的相互作用示意。绿色为holo/active晶体结构，橙色为holo/inactive初始模型。A展示“capping triad”内的盐桥网络，E展示Y480–R528氢键。 (B, F) 对应A与E的距离分布箱线图，绿色为holo/active N3Aout/BCin集合，橙色为holo/inactive N3Ain/BCout集合，绿色/红色线标记晶体结构与初始模型的距离。绿色箱体（左）表示接触更短更稳，橙色（右）表示接触被拉开。 (C, D) 来自MD轨迹的代表性结构，进一步对比“capping triad”的几何组合。active集合保持三联体稳定相互作用，而inactive集合中Q532更倾向远离R484，仅保留D533与R484的单盐桥。相互作用类型 Holo/Active状态 Holo/Inactive状态结构后果 R484-Q532盐桥稳定存在（绿色箱体分布靠左）被破坏/不稳定（橙色箱体分布右移） Q532远离R484，triad结构解体 R484-D533盐桥稳定存在相对保持（单盐桥） D533靠近R484，但Q532已远离 Y480-R528氢键稳定存在显著减弱 αC螺旋与PBC的空间解耦这些差异与文献中的突变结果一致，支持用holo/active与holo/inactive两组MD集合来筛选激活所必需的PBC/αC螺旋接触。因此在N3Ain/BCout集合中，这些接触应被明显削弱，而在N3Aout/BCin集合中保持稳定，这正是B–F所观测到的趋势。 (G–J) R484A突变体的CHESPA分析：G为矢量示意，H为WT与R484A在cGMP结合状态下的化学位移差异，I为fractional shift（$X$），J为$\cos(\Theta)$。CHESPA用WT的apo→holo位移变化定义激活向量，用突变体相对WT的位移变化定义突变向量，比较方向与投影大小。激活向量由WT在apo与holo之间的化学位移差值组成，代表配体结合引发的构象变化方向。这些化学位移来自实验NMR 1H–15N HSQC谱图，在WT与R484A的apo与cGMP结合条件下测量后进行CHESPA投影分析。 $\cos(\Theta)$计算式： \[\cos(\Theta)=\frac{\vec{\delta}_{\text{mut}}\cdot\vec{\delta}_{\text{act}}}{\left|\vec{\delta}_{\text{mut}}\right|\left|\vec{\delta}_{\text{act}}\right|}\] $X$值计算式： \[X=\frac{\vec{\delta}_{\text{mut}}\cdot\vec{\delta}_{\text{act}}}{\left|\vec{\delta}_{\text{act}}\right|^{2}}\] $X$表示突变效应在激活方向上的投影强度，$X=0$表示不沿激活方向变化，$X<0$说明突变把体系拉回非活化方向。 Δδ表示综合化学位移差异强度，用于衡量突变对局部结构的总体扰动幅度。多数残基$X$为负且$\cos(\Theta)$接近−1，说明R484A显著把体系拉回非活化方向，验证R484是维持active构象的关键锚点。 Capping triad是PfPKG CBD-D激活的关键结构元件，由PBC的R484与C端αC螺旋的Q532/D533形成的盐桥网络组成。这一结构在PfPKG中是独特的，人类PKG Iβ对应位置是K308，不与αC螺旋形成类似的相互作用（Figure S5），这为设计物种选择性抑制剂提供了基础。 R484的位置优势：R484位于PBC loop，其guanidinium基团可以同时与Q532和D533形成离子对立体化学互补：在active构象中（N3Aout/BCin），R484、Q532、D533三者空间排列形成稳定的三角网络双重稳定作用：Capping triad既稳定了αC螺旋的向内构象（BCin），又通过R484-cGMP cation-π相互作用稳定了配体结合 4. cGMP结合相互作用：激活与非活性态的差异进一步分析cGMP与PBC和BBR区域的相互作用，可以识别区分holo/active和holo/inactive状态的关键结合特征。图7：PBC与cGMP及类似物的关键相互作用 (A–C) cGMP与PfPKG CBD-D的相互作用示意（PDB: 4OFG），虚线标示监测的相互作用距离，标注参与相互作用的残基 (D, E) 关键原子对距离分布的箱线图，绿色为holo/active N3Aout/BCin，橙色为holo/inactive N3Ain/BCout，红色虚线框标示两种集合间变化最显著的相互作用 (F–H) 磷酸硫代cGMP类似物的结构示意：Sp-cGMPS和Rp-cGMPS (I) PfPKG 401-853的环核苷酸依赖性激活曲线，展示不同类似物的激活能力 Figure 7A-C详细展示了cGMP如何与PBC和BBR区域形成多重相互作用：区域 cGMP部分关键残基相互作用类型功能 PBC 磷酸基团 482-485, 492-493 氢键网络锚定cGMP的磷酸基团 PBC 磷酸基团 T493 桥接氢键连接轴向氧和氨基 BBR 鸟嘌呤碱基 R473 氢键识别碱基特异性 PBC 鸟嘌呤碱基 R484 cation-π 稳定碱基结合，形成capping triad的一部分 T493的羟基同时与cGMP的磷酸基团（轴向氧）和氨基形成氢键，在空间上起到桥梁作用，是PBC区域中唯一同时与cGMP两个部分相互作用的残基。Figure 7D, E的红色虚线框标出了两种holo状态间差异最大的相互作用： A485-cGMP氢键：Holo/active中稳定，holo/inactive中被破坏（Figure 7D） R484-cGMP cation-π相互作用：Holo/active中强，holo/inactive中显著减弱（Figure 7E）这两个相互作用的变化与Figure 6中Capping triad的破坏相呼应，共同导致了holo/inactive状态的失活。 cGMP类似物的设计策略与实验验证文献基于MD预测设计了Rp-cGMPS和Sp-cGMPS两种立体异构体，用于验证A485-cGMP氢键的重要性：类似物修饰位置设计原理预测效果实验结果 Rp-cGMPS (Figure 7H) 轴向氧→硫（Rp构型）破坏A485-cGMP关键氢键激酶活性大幅降低 75%活性降低，验证预测 Sp-cGMPS (Figure 7G) 平分向氧→硫（Sp构型）修饰非关键相互作用活性轻微降低仅10%降低，作为对照 Figure 7I的激酶活性实验显示，Rp-cGMPS的弱激动剂效应（蓝色曲线）激活能力降至~25%，证明A485-cGMP氢键对激酶激活至关重要；Sp-cGMPS的部分激动剂效应（黑色曲线）激活能力降至~90%，验证了其他相互作用的保守性。这形成了从预测到验证的闭环：MD模拟（Figure 7D, E）→设计类似物→激酶活性实验（Figure 7I）。变构抑制剂的启示 Figure 7的结果揭示了靶向R484-A485-cGMP相互作用网络的潜力：选择性破坏：这两个相互作用在holo/active中强，在holo/inactive中弱，是理想的变构抑制剂靶点保留结合亲和力：其他cGMP-PBC/BBR相互作用在两种holo状态中保守，破坏R484-A485不会完全丧失cGMP结合物种选择性基础：PfPKG的R484可形成capping triad，而人类PKG Iβ的K308不与αC螺旋相互作用（Figure S5），为宿主-寄生虫选择性提供了结构基础唉，其实这些都是如何解释机制能算的一些指标。虽然都能用，但是似乎还是没有那么直接，比如直接去算QM过程的free energy vs RC。讨论本研究通过MD模拟完整映射了PfPKG CBD-D的四态变构热力学循环，识别了区分激活与非活性状态的关键相互作用。这些发现为理解PfPKG的变构调控机制提供了原子级视角。变构抑制剂设计的结构基础 holo/inactive中间态代表了配体结合但不激活的独特状态，是设计变构抑制剂的关键靶点。通过比较holo/active和holo/inactive状态，我们识别了几个关键的结构差异：关键相互作用 Holo/Active状态 Holo/Inactive状态抑制剂设计策略 R484-cGMP阳离子-π作用强（稳定）弱或缺失设计类似物削弱此作用 A485-cGMP氢键完整（氧原子）破坏 Rp-cGMPS中氧→硫替代显著降低活性 R484-Q532/D533-capping triad 存在弱化或缺失靶向破坏此三联体 C端螺旋-αC螺旋相互作用稳定松动设计分子阻止螺旋靠近 Rp-cGMPS的实验验证将A485酰胺与cGMP磷酸氧的氢键破坏后（氧→硫替代），激酶活性降低75%，证明了靶向R484-A485相互作用可以实现变构抑制，同时保持与cGMP其他接触的保守性。物种选择性机制 PfPKG的R484可形成capping triad与C端αC螺旋的Q532/D533相互作用，而人类PKG Iβ对应的K308不与αC螺旋相互作用（Figure S5）。靶向R484相互作用可能实现PfPKG vs人类宿主的选择性。 Q&A Q1：为什么PBC区域的稳定化需要同时满足cGMP结合和变构构象变化？ A1：PBC区域的动力学响应显示出独特的双重依赖机制，这在物理化学上可以通过以下几个方面理解：构象选择的局限性：如果纯粹是构象选择机制（蛋白预先存在multiple conformations，cGMP选择其中一种结合），那么apo/active状态（已经具有active构象）的PBC应该也相对稳定。但Figure 3B和4B显示，apo/active的PBC RMSF和RMSD都显著高于holo/active，说明仅有active构象是不够的。诱导契合的局限性：如果纯粹是诱导契合机制（cGMP结合后诱导蛋白构象改变），那么holo/inactive状态（有cGMP结合）的PBC应该相对稳定。但数据显示holo/inactive的PBC RMSF和RMSD与apo/inactive相近，说明仅有cGMP结合也是不够的。协同作用的物理本质：cGMP与PBC的相互作用形成一个正反馈循环： cGMP优先结合到active构象的PBC（构象选择成分）：active构象的PBC具有更适合的几何形状和电荷分布，结合亲和力更高 cGMP结合进一步稳定和锁定active构象（诱导契合成分）：cGMP与PBC的氢键、cation-π等相互作用网络增强了active构象的稳定性这两个过程是同时发生、相互促进的，而非先后独立的步骤能量景观的视角：在四态热力学循环中，holo/active状态位于能量最低点（Figure 5的右上象限聚集了大量数据点），而apo/active和holo/inactive都位于较高的能量状态。这表明cGMP结合和active构象的同时满足才能达到最稳定的能量状态，两者存在协同的能量贡献。 Q2：为什么所有激活路径都必须经过“紧凑”N3Ain/BCin中间态？ A2：这一发现可以通过能量景观理论和拓扑约束来解释：拓扑约束的物理原因：从N3Ain/BCout（inactive）到N3Aout/BCin（active）的转变涉及两个主要结构元件的重排。直接从N3Ain/BCout跳变到N3Aout/BCin需要同时改变N3A和αB-αC的位置，这在能量上是不利的。相反，通过紧凑的N3Ain/BCin中间态，可以逐步改变各个元件的位置，降低能垒。 N3A的in/out采样动力学：Figure 5显示N3A在所有四个状态下都能动态采样in和out取向，这意味着N3A的重排相对容易。而αB-αC螺旋的in/out转变则更受构象状态的约束（Figure 4C显示active状态αB-αC更稳定）。因此，N3Ain/BCin中间态代表了一个能量上的有利过渡态，其中N3A已经向内，αB-αC也准备向内移动。与HCN通道的比较：HCN通道的变构转变遵循多分支路径，而PfPKG CBD-D显示出对紧凑中间态的强偏好，这反映了不同环核苷酸结合结构域的变构调控机制多样性，可能与功能需求（如激活速度、调控精度）相关。 Q3：holo/inactive中间态如何指导变构抑制剂设计？ A3：holo/inactive中间态代表了配体结合但不激活的独特状态，其结构特征为设计变构抑制剂提供了三个关键策略：靶向R484-A485与cGMP相互作用：Figure 7D, E显示从holo/active到holo/inactive转变时，R484-cGMP的cation-π相互作用和A485-cGMP氢键被显著破坏。Rp-cGMPS实验（Figure 7I）证明破坏A485-cGMP氢键可降低75%激酶活性，这验证了靶向这些相互作用可以实现变构抑制。破坏capping triad相互作用：Figure 6显示R484与C端αC螺旋的Q532/D533形成的capping triad在holo/active状态稳定存在，而在holo/inactive状态被破坏。设计小分子或肽段干扰这个三联体，可以阻止C端螺旋与PBC的稳定相互作用，从而抑制激活。物种选择性的结构基础：Figure S5显示PfPKG的R484可形成capping triad与C端αC螺旋相互作用，而人类PKG Iβ对应的K308不与αC螺旋形成类似相互作用。这种差异为设计PfPKG选择性抑制剂提供了明确靶点，可以实现对疟原虫的选择性毒性，避免对人类宿主的副作用。关键结论与批判性总结主要结论本研究的结论与原文讨论部分一致，可归纳为以下几点：完整描绘四态热力学循环的动力学变化：通过MD与实验数据支持，系统刻画了apo/inactive、apo/active、holo/inactive、holo/active四态的动力学差异，尤其涵盖实验难以直接表征的中间态。区分cGMP结合与变构构象变化的贡献：动力学地图揭示apo/inactive→holo/active转变同时依赖cGMP结合与构象变换，两者贡献可被拆分比较。中间态的结构特征具有设计价值：相似性分析显示apo/active兼具apo/inactive与holo/active特征，holo/inactive更接近apo/inactive，这为“结合但不激活”的变构抑制剂提供了明确参照。关键接触位点明确：PBC与αC螺旋的接触（R484‑Q532/D533 capping triad、Y480‑R528氢键）对激活至关重要，且R484‑A485与cGMP的相互作用在holo/inactive与holo/active之间差异显著，提示可优先靶向这些接触进行选择性干预。物种选择性线索：PfPKG的R484对应人类PKG Iβ的K308，后者不与αC螺旋形成同类接触，破坏R484相关相互作用可能带来Pf与宿主的选择性。已知限制与待验证点中间态的实验表征仍具挑战：原文指出apo/active与holo/inactive等中间态难以通过实验直接捕捉，因此目前主要依赖模拟与间接实验证据支撑。研究意义与可预期方向变构抑制剂设计的直接线索：holo/inactive特征可用于设计“结合但不激活”的配体，优先削弱R484‑A485与cGMP的作用或破坏capping triad。验证路径清晰：文中通过突变与CHESPA证实R484A可逆转激活方向，支持以PBC/αC螺旋接触为核心的验证与优化策略。
Molecular Dynamics · 2026-01-22
LSP-MD：捕捉热振动驱动变构效应的快速计算方法
 LSP-MD：捕捉热振动驱动变构效应的快速计算方法本文信息标题：LSP-MD： A Fast Computational Method to Study Allostery Driven by Thermal Vibrations 作者：Alexandr P. Kornev 发表时间： 2025年11月4日单位：LSP Consulting LLC（美国加利福尼亚州）引用格式：Kornev, A. P. (2025). LSP-MD: A Fast Computational Method to Study Allostery Driven by Thermal Vibrations. Journal of Chemical Theory and Computation, 21(21), 8699-8710. https://doi.org/10.1021/acs.jctc.5c01094 源代码/软件：论文未公开代码，但LSP Consulting LLC提供与LSP相关方法的咨询服务和许可证（见Conflict of Interest声明）摘要与热振动相关的构象熵在蛋白质功能中发挥根本性作用，从配体结合和催化到变构调节。Cooper和Dryden首次将熵驱动变构作为这些效应的一个例子提出。然而，测量底层热运动在技术上仍然具有挑战性。在此，我们介绍了LSP-MD，这是一种建立在局部空间模式（LSP）对齐基础上的计算方法，用于跟踪分子动力学（MD）模拟中的侧链稳定性。LSP-MD使用基于图的蛋白质残基网络（PRNs），其边权重来源于快速的局部几何涨落。应用于蛋白激酶A（PKA）时，该方法捕获了皮秒时间尺度的振动，振幅在0-2Å范围内，波数低于100 $\mathrm{cm^{-1}}$，正好在熵介导信号传导的范围内。从LSP-MD网络导出的中心性指标在不同模拟长度、向量定义和力场下保持稳定，确认了鲁棒性。重要的是，LSP-MD重现了传统LSP分析的关键发现，同时提供了更清晰的物理基础和更高的计算效率。该方法为探索各种大分子系统中的熵驱动变构行为开辟了新机会。核心结论热振动的直接测量：LSP-MD方法首次实现了对皮秒时间尺度热振动的直接量化，捕获了振幅0-2Å、波数低于100 $\mathrm{cm^{-1}}$的振动模式网络化稳定性分析：通过基于蛋白质残基网络（PRN）的中心性指标，将局部几何涨落转化为全局变构信号计算效率提升：相比传统LSP对齐方法，LSP-MD消除了耗时的模式搜索和结构映射步骤，可将500帧轨迹分析，而传统方法仅能处理100帧方法鲁棒性验证：中心性指标在不同模拟长度（10-100 ns）、采样率、向量定义和力场（ff14SB与CHARMM36）下保持高度稳定物理意义明确：用单一物理参数Z（几何偏差的欧几里得范数）量化残基对稳定性，替代了传统方法的ad hoc参数背景蛋白质在沿着折叠漏斗向其天然结构滑动时，随着结构变得更加有序，其熵会减少。然而，即使在折叠完成后，侧链仍然保留了相当大的流动性。这种残留熵，也称为构象熵，在蛋白质功能中发挥着重要作用。在他们最近的综合综述中，Wankowicz和Fraser证明这些熵效应是蛋白质动力学的普遍特征，影响着从配体结合特异性到酶催化、从蛋白质稳定性到变构信号传导的各个方面。这些效应在变构调节中尤其重要，其中配体在一个位点的结合会通过结构变化或动力学效应远程影响另一个位点的功能。早在1984年，Cooper和Dryden就提出了一个革命性的概念：蛋白质的变构效应可以完全由熵变化驱动，而不需要明显的结构重排。他们计算表明，侧链构象熵的微小变化（每个残基约0.4-1.2 kJ/mol）就足以产生显著的变构效应。这一预测在过去几十年中得到了实验支持。核磁共振（NMR）弛豫测量、异核核Overhauser效应和顺序参数分析等实验技术已经能够直接探测这些快速的热运动。然而，这些实验方法通常需要昂贵的设备、专业的样品制备（如同位素标记），并且难以获得全原子级别的分辨率。从计算角度看，分子动力学（MD）模拟提供了研究这些热振动的理想工具。现代MD模拟可以在飞秒时间分辨率下跟踪每个原子的运动，理论上可以捕获从皮秒到毫秒时间尺度的所有动力学过程。然而，从海量轨迹数据中提取有意义的变构信号仍然是一个巨大的挑战。传统的分析方法要么过于简化（如均方根偏差分析），要么计算成本过高（如全原子互相关分析）。为了解决这个问题，Kornev等人此前开发了局部空间模式（LSP）对齐方法，用于比较蛋白质晶体结构并识别侧链稳定性的变化。LSP方法通过将残基表示为向量，并分析不同结构中残基对之间几何关系的变化，成功捕获了与变构相关的稳定性模式。然而，传统LSP方法依赖于大量晶体结构的比较，且需要进行穷举式的模式搜索和结构映射，计算成本高昂，限制了其在MD轨迹分析中的应用。关键科学问题热振动的量化难题：如何从MD模拟的海量轨迹数据中提取出真正与变构相关的微小热振动信号，而不是被其他大尺度构象变化所淹没时间尺度的匹配问题：变构相关的热振动主要发生在皮秒到纳秒时间尺度，如何设计专门针对这一时间尺度的高效分析方法物理意义的阐释：如何将抽象的网络拓扑参数与具体的物理过程（热振动、构象熵）联系起来，提供明确的物理解释计算效率与准确性的平衡：如何在保持对变构信号敏感的同时，大幅降低计算成本，使方法能够应用于大规模的MD轨迹分析创新点 LSP-MD方法框架：提出了一种全新的MD轨迹分析方法，直接在轨迹内量化残基对的稳定性，无需与外部参考结构比对 Z参数的引入：使用几何偏差的欧几里得范数作为单一稳定性指标，具有明确的物理意义，替代了传统LSP方法的ad hoc参数网络化变构分析：将局部稳定性信息转化为PRN的边权重，通过网络中心性指标（DC、BC）识别关键的变构节点系统性的参数优化：系统研究了模拟时间、样本大小、距离截断等参数对结果的影响，提供了标准化的分析流程方法验证与对比：与传统LSP对齐方法进行了系统对比，证明新方法不仅计算效率更高，而且保留了原有的核心发现研究内容 LSP-MD方法的原理与实现！fig1 图1：LSP-MD方法的局部稳定性测量原理该图展示了LSP-MD如何通过四个几何距离量化残基对稳定性： (A) 蛋白质残基网络（PRN）示意图，节点为残基，边的粗细反映稳定性权重 (B) 残基向量化几何定义，展示两个残基向量间的四个距离（$d_1, d_2, d_3, d_4$） (C) Z参数计算流程：四个距离偏差（$\Delta d_1, \Delta d_2, \Delta d_3, \Delta d_4$）通过欧几里得范数组合为Z (D) PKA系统的距离偏差分布散点图，蓝色点为标准向量，红色点为长侧链向量，展示Z值集中在0-2 Å范围 Scheme 1：LSP对齐方法与LSP-MD算法的流程对比该图对比了传统LSP对齐方法和LSP-MD方法的计算流程： (A) LSP对齐算法：用于比较两个不同的蛋白质结构。首先计算两个蛋白质中所有残基对的内部几何关系，然后进行计算密集型的相似性搜索（红色矩形标注），寻找两个蛋白质中具有相似空间模式的残基对。最终输出一组同构子图，显示两个蛋白质中的相似模式 (B) LSP-MD算法：用于分析单个蛋白质在多个构象下的动力学特征。对轨迹中的每一帧计算所有残基对的内部几何关系，然后对整个轨迹取平均，计算几何偏差，最终得到稳定性指标（Z值）。输出单一的PRN图，表征蛋白质的构象动力学关键区别：传统LSP需要在两个蛋白质之间进行穷举式的模式搜索（计算复杂度高），而LSP-MD只需在单个蛋白质的轨迹内计算平均和偏差（计算效率高）。LSP-MD用时间平均替代了结构比对，用几何涨落替代了模式相似性。核心思想：从几何涨落到网络权重 LSP-MD的核心思想是将MD轨迹中每个残基对的局部几何稳定性量化为一个单一的物理参数，然后将其转化为蛋白质残基网络（PRN）的边权重，通过网络分析识别关键的变构节点。方法的具体实现步骤 1。残基向量化：将每个残基表示为一个向量，通常从Cα指向Cβ。对于甘氨酸（没有Cβ）或其他特殊情况，可以使用替代定义（如N-Cα或质心-Cα） 2。距离定义：对于两个残基的向量对（残基 $i$ 的向量为$\mathbf{v}_i$，残基 $j$ 的向量为$\mathbf{v}_j$），定义四个距离： $d_1$：残基 $i$ 的起点到残基 $j$ 的起点 $d_2$：残基 $i$ 的起点到残基 $j$ 的终点 $d_3$：残基 $i$ 的终点到残基 $j$ 的起点 $d_4$：残基 $i$ 的终点到残基 $j$ 的终点 3。轨迹平均：计算整个MD轨迹中这四个距离的平均值$\langle d_1 \rangle, \langle d_2 \rangle, \langle d_3 \rangle, \langle d_4 \rangle$ 几何偏差计算：对于轨迹中的每一帧，计算四个距离的偏差$\Delta d_k = d_k - \langle d_k \rangle$（$k=1,2,3,4$） Z参数计算：将四个偏差组合为单一参数Z，使用欧几里得范数： \(Z = \sqrt{(\Delta d_1)^2 + (\Delta d_2)^2 + (\Delta d_3)^2 + (\Delta d_4)^2}\) 边权重转换：将Z值转换为边权重W，使用公式$W = \exp(-Z)$。这样稳定的残基对（小Z）获得高权重，不稳定的残基对（大Z）获得低权重网络构建：仅当两个残基的Cα原子距离小于截断值（通常为12Å）时，在它们之间创建边中心性分析：计算加权PRN中每个节点的度中心性（DC）和介数中心性（BC），识别关键的变构节点 graph TB Start["MD轨迹输入"] --> S1 subgraph S1["1.残基向量化"] direction LR A1["Cα→Cβ向量定义"] --> A2["替代向量定义<br/>甘氨酸/末端残基"] end S1 --> S2 subgraph S2["2.几何参数提取"] direction LR B1["定义4个距离<br/>d1, d2, d3, d4"] --> B2["计算轨迹平均<br/>⟨d⟩值"] --> B3["计算偏差<br/>Δd = d - ⟨d⟩"] end S2 --> S3 subgraph S3["3.稳定性量化"] direction LR C1["计算Z参数<br/>欧几里得范数"] --> C2["转换为边权重<br/>W = exp(-Z)"] end S3 --> S4 subgraph S4["4.网络构建与分析"] direction LR D1["构建PRN<br/>Cα距离<12Å"] --> D2["计算DC和BC<br/>识别关键节点"] end S4 --> Result["输出变构热点图谱"] Z参数的物理意义 Z参数是LSP-MD方法的核心创新，它具有明确的物理意义：几何稳定性的直接度量：Z值反映了残基对之间相对几何关系偏离其轨迹平均状态的程度。小Z值表示残基对的相对位置保持稳定，大Z值表示几何关系波动较大热振动幅度的表征：在PKA的10纳秒模拟中，Z值主要分布在0-2Å范围内，这与热振动引起的小幅度构象变化一致波数选择性：通过快速傅里叶变换（FFT）分析发现，Z值变化的波数分量主要集中在100 $\mathrm{cm^{-1}}$以下，正好对应于热激发模式的波数范围（<200 $\mathrm{cm^{-1}}$）与传统LSP对齐方法的区别传统LSP对齐方法需要比较多个实验结构（通常是不同配体结合状态的晶体结构），通过穷举式的模式搜索和结构映射来识别侧链稳定性的变化。LSP-MD方法与传统LSP方法的关键区别总结如下：特征传统LSP对齐方法 LSP-MD方法数据来源需要多个高质量晶体结构（不同配体状态）直接在MD轨迹内分析，无需外部参考结构计算成本模式搜索和结构映射耗时长，难以处理大量轨迹消除模式搜索和结构映射，计算效率显著提升参数设置使用ad hoc阈值参数，物理意义不明确使用Z参数（几何偏差的欧几里得范数），物理意义明确适用范围受限于可获得晶体结构的系统可应用于任何MD模拟系统处理规模通常限于100帧左右结构对比可轻松处理500帧甚至更多轨迹帧应用案例：蛋白激酶A的热振动分析系统选择与模拟设置蛋白激酶A（PKA）是研究变构调节的经典模型系统。PKA具有典型的双叶激酶折叠，包括较小的N叶（主要包含β折叠）和较大的C叶（主要包含α螺旋）。两叶之间的铰链区域包含了催化位点和多个关键的调节元件，如glycine-rich loop和αC-螺旋。研究者使用PKA的催化亚基进行测试，模拟设置总结如下：参数类别具体设置说明/目的初始结构 PDB ID 1ATP ATP结合状态的PKA催化亚基力场 AMBER ff14SB 蛋白质标准力场溶剂模型 TIP3P水，10Å缓冲水化蛋白，提供真实溶剂环境离子条件 Na⁺/Cl⁻，150 mM 中和电荷，模拟生理盐浓度平衡协议逐步加热至300 K，1 atm 系统平衡至目标温度和压强生产模拟 10 ps（0.5 fs步长）高分辨率轨迹，捕获皮秒振动 10-100 ns（2 fs步长）常规轨迹，稳定性分析模拟软件 AMBER 20 皮秒时间尺度的热振动特征图2：PKA中代表性残基对的Z值时间演化与频谱分析该图从多个时间尺度展示了LSP-MD捕获的热振动特征： (A) 皮秒时间尺度的Z值演化（1 ps轨迹，0.5 fs步长）：曲线展示了三个代表性残基对的Z值随时间的超精细变化。黑色曲线（K72-E91）：连接N叶β折叠和调节性αC-螺旋的保守盐桥，被视为激酶活性态的标志。曲线非常平滑，Z值变化极小（千分之一埃量级），展现了极高的结构刚性红色曲线（I150-D220）：位于C叶内部的残基对，Z值略高于盐桥，反映了相对温和的灵活性蓝色曲线（G55-G186）：连接glycine-rich loop和DFG基序的残基对，Z值变化最为明显，代表了分子中最可动的区域插图：三个残基对在PKA结构上的位置。较大的C端用棕褐色着色，清晰显示了两叶结构和铰链区域这些超精细轨迹显示了LSP-MD方法的时间分辨率优势：即使在0.5 fs步长下，Z值曲线仍然非常平滑，能够捕捉到残基运动的每一个细节。 (B) K72-E91盐桥Z值变化的频谱分析：通过快速傅里叶变换（FFT）将时域信号转换为频域功率谱。横轴为波数（$\mathrm{cm^{-1}}$），纵轴为相对功率（%）。关键发现：主波数分量集中在100 $\mathrm{cm^{-1}}$以下，最高功率谱峰出现在6.6 $\mathrm{cm^{-1}}$（>12%相对功率）。这一低频分布正好对应于热激发模式的波数范围（<200 $\mathrm{cm^{-1}}$），证明了LSP-MD捕获的振动确实是由热运动驱动的。这一波数分布具有双重意义：低于热激发阈值：蛋白质中可以热激发的振动模式波数阈值约为200 $\mathrm{cm^{-1}}$。LSP-MD捕获的振动波数（5-100 $\mathrm{cm^{-1}}$）完全在这一范围内，说明这些振动确实是由热运动驱动的与变构相关的波数范围：先前研究表明，小的变构事件（如侧链重新取向）主要影响100 $\mathrm{cm^{-1}}$以下的低波数模式。LSP-MD正是聚焦于这一关键的波数窗口 (C) 纳秒时间尺度的Z值演化（100 ns轨迹）：展示了更长时间尺度下Z值的变化。蓝色曲线（G55-G186）：Z值最大可达约5Å，出现多个峰，对应于glycine-rich loop的大幅度构象重排黑色和红色曲线（K72-E91和I150-D220）：Z值变化相对温和，最大约3Å，反映了刚性结构域的稳定性视觉检查发现，这些Z值的峰值对应于构象状态的转变，如loop的闭合/开放、侧链的rotameric跳跃等。 (D) 不同长度模拟的Z值分布统计：直方图展示了从不同长度模拟（100 ps、1 ns、10 ns、100 ns）中提取的500个PKA结构中所有残基对的Z值频率分布。横轴为Z值（Å），右端点表示Z>2Å的统计。 10 ns模拟：Z值主要集中在0-1Å范围 100 ns模拟：分布略微变宽，但绝大多数残基对的Z值仍低于2Å 这一发现表明，尽管存在可动区域（如loop），PKA的大部分残基对在纳秒时间尺度上仍然保持着相对稳定的几何关系。这种局部稳定性是蛋白质三维结构得以维持的基础，也是LSP-MD方法能够捕获有意义信号的前提。模拟时间对中心性指标的影响研究者系统地研究了模拟时间对度中心性（DC）和介数中心性（BC）的影响：图3：模拟时间对LSP-MD中心性指标的影响该图系统展示了不同模拟长度下LSP-MD网络的收敛行为： (A) 度中心性（DC）随模拟时间的变化：折线图展示了αF-螺旋中12个连续残基的DC值在不同模拟长度下的变化（误差棒为5次独立重复的标准误差）。关键发现：在10 ns之前，DC值明显被高估，随后快速下降并趋于平稳。这表明短暂模拟（<10 ns）未能充分探索热振动的完整范围，导致边权重整体偏高 (B) 介数中心性（BC）随模拟时间的变化：同样的12个αF-螺旋残基的BC值变化。关键发现：与DC相反，BC值在短模拟中被低估，随模拟时间增加而上升。这是因为BC对全局网络拓扑更敏感，短模拟中的高边权重掩盖了真实的通信路径结构 (C) 所有残基DC值的标准误差分布：箱线图展示了PKA全部338个残基在不同模拟时间下DC值的重复性（5次重复的标准误差）。横轴为模拟长度，纵轴为标准误差。关键发现：标准误差在达到10 ns后基本稳定，更长的模拟并不会显著增加噪声 (D) 所有残基BC值的标准误差分布：与DC类似，BC的标准误差也在10 ns后收敛。注意：BC的绝对误差值高于DC，这与BC对全局网络结构的敏感性一致 (E) 10 ns与100 ns模拟的DC值相关性：散点图对比了所有残基在这两种模拟长度下的DC值。Pearson相关系数$r=0.997$，表明极高的一致性。大多数点沿对角线紧密分布，说明10 ns和100 ns的DC图谱几乎相同 (F) 10 ns与100 ns模拟的BC值相关性：BC值的对比也显示出强相关性（$r=0.987$），虽然略低于DC，但仍证明10 ns模拟已能捕获关键的变构通信路径中心性指标的定义在详细讨论结果之前，我们先明确两个核心网络分析指标的定义和物理意义：度中心性（Degree Centrality, DC）衡量节点在网络中的直接连接重要性。在加权PRN中，节点 $i$ 的DC定义为与该节点相连的所有边的权重之和： \[\mathrm{DC}(i) = \sum_{j \in N(i)} W_{ij}\] 其中 $N(i)$ 是节点 $i$ 的邻居集合，$W_{ij} = \exp(-Z_{ij})$ 是节点 $i$ 和 $j$ 之间的边权重。DC反映了一个残基与周围残基形成稳定连接的能力。高DC残基通常位于蛋白质结构的稳定核心，与其周围的残基保持紧密且稳定的几何关系。介数中心性（Betweenness Centrality, BC）衡量节点在网络中作为”桥梁”或”中继”的能力。节点 $i$ 的BC定义为： \[\mathrm{BC}(i) = \sum_{s \neq i \neq t} \frac{\sigma_{st}(i)}{\sigma_{st}}\] 其中 $\sigma_{st}$ 是从节点 $s$ 到节点 $t$ 的最短路径总数，$\sigma_{st}(i)$ 是经过节点 $i$ 的最短路径数。BC反映了残基在网络通信中的重要性。高BC残基通常位于不同结构域之间的通信路径上，充当变构信号的”中继站”，在长距离信号传导中发挥关键作用。这两个指标共同刻画了残基在蛋白质变构网络中的角色：DC反映局部稳定性，BC反映全局通信能力。 10 ns模拟时间转折点分析模拟时间 DC值表现 BC值表现收敛状态物理原因 <10 ns 被高估被低估未收敛未能充分探索热振动范围，$\langle d \rangle$偏向起始构象，导致$\Delta d$偏小，Z值偏低，边权重偏高 ≥10 ns 趋于稳定趋于稳定充分收敛 $\langle d \rangle$已充分收敛，DC和BC标准误差稳定，10 ns与100 ns相关性$r>0.98$ 这一发现的实际意义是：对于PKA这类蛋白质，10 ns模拟已足够捕获热振动驱动的变构信号，更长的模拟并不会显著改变中心性图谱。这大大降低了计算成本，使LSP-MD方法能够应用于大规模的蛋白质动力学研究。样本大小的优化除了模拟时间，研究者还研究了从轨迹中采样的帧数对结果的影响：图4：样本大小对LSP-MD中心性指标的影响该图评估了从10 ns轨迹中提取不同数量帧对分析结果的影响： (A) DC值随样本大小的变化：折线图展示了αF-螺旋中12个残基的DC值随采样帧数增加的变化（从5帧到2500帧）。横轴为帧数（对数坐标），纵轴为DC值。关键发现：DC值在小样本（<100帧）时波动较大，在约100帧时趋于稳定 (B) 所有残基DC值的标准误差分布：箱线图展示了PKA全部338个残基在不同样本大小下DC值的重复性（5次重复的标准误差）。关键发现：标准误差随样本增加而下降，在约100-500帧时达到平台期 (C) BC值随样本大小的变化：同样的12个αF-螺旋残基的BC值变化。BC值需要更多帧才能收敛，反映了其对全局网络结构的敏感性 (D) 所有残基BC值的标准误差分布：BC的标准误差在约500帧时达到较好的稳定性 (E) 100帧与2500帧的DC值相关性：散点图对比了这两种采样密度的DC值。Pearson相关系数$r=0.98$，说明100帧已能代表完整轨迹的DC图谱 (F) 100帧与2500帧的BC值相关性：BC值的相关性（$r=0.96$）同样很高，证明约100帧的采样已足够使用10 ns轨迹（每4 ps保存一帧，共2500帧），不同采样帧数的性能对比：采样帧数 DC和BC稳定性计算开销推荐程度 <100帧波动较大，标准误差高低不推荐 ~100帧趋于稳定低可接受 500帧提供更好的稳定性小推荐建议的平衡方案是使用约500帧进行分析。考虑到LSP-MD的高效性，处理500帧的计算时间非常短，这一建议具有很高的实用性。距离截断的优化 PRN的构建需要定义一个距离截断，只有两个残基的Cα原子距离小于该截断值时才创建边。研究者系统测试了不同截断值的影响：图5：Cα距离截断对LSP-MD网络拓扑的影响该图系统评估了不同距离截断值对PRN结构和中心性指标的影响： (A) 不同截断距离下的ForceAtlas2网络布局：使用力导向算法可视化PRN拓扑结构，节点大小反映DC，颜色深浅反映BC。展示了从8Å到16Å截断的网络密度和模块化程度变化 (B) 模块化和边密度随截断距离的变化曲线：绿色曲线（模块化）：衡量网络划分为内部凝聚模块的能力。纵轴为模块化指数，横轴为截断距离。关键发现：在10-15Å范围出现明显的斜率变化（红色虚线标注），二阶差分（插图）确认了12Å是最优截断值蓝色曲线（边密度）：实际边数与可能的最大边数之比。边密度随截断增加而单调上升，但在10-15Å范围出现斜率变化 (C) 12Å与14Å截断的DC值相关性：散点图对比了这两种截断下所有残基的DC值。Pearson相关系数$r=0.96$，说明在12-14Å范围内DC值高度一致，网络拓扑保持稳定 (D) 12Å与14Å截断的BC值相关性：BC值的相关性（$r=0.86$）同样显著，证明了这一截断范围的鲁棒性网络拓扑的变化截断距离网络特征模块化程度连通性适用性 8 Å 网络非常稀疏，节点分散高差不推荐 10 Å 网络开始形成基本骨架较高较差可接受 12 Å 网络密度适中，模块清晰可见，高BC节点集中在模块中心稳定良好推荐 14 Å 网络进一步致密化，模块边界开始模糊适中很好可接受 16 Å 网络非常密集显著下降过度连通不推荐定量指标含义模块化指数（Modularity Q）衡量网络划分为内部凝聚模块的程度，定义为： \(Q = \frac{1}{2m} \sum_{i,j} \left[ W_{ij} - \gamma \frac{k_i k_j}{2m} \right] \delta(c_i, c_j)\) 其中： $W_{ij}$ 是节点 $i$ 和 $j$ 之间的边权重（在LSP-MD中为 $\exp(-Z_{ij})$） $k_i = \sum_j W_{ij}$ 是节点 $i$ 的加权度 $m = \frac{1}{2} \sum_{i,j} W_{ij}$ 是网络中所有边的权重总和 $\gamma$ 是分辨率参数（通常为1） $\delta(c_i, c_j) = 1$ 如果节点 $i$ 和 $j$ 在同一模块，否则为0 如何理解模块化指数？用一个社交网络类比：模块化指数Q衡量网络能否清晰地分成几个内部紧密、外部疏离的“小圈子”。计算逻辑（简化版）： \(Q \approx \frac{\text{圈子内部的实际联系数} - \text{随机期望的内部联系数}}{\text{总联系数}}\) Q接近1（高度模块化）：三个完全不交流的微信群（科研群、游戏群、购物群），群内互动频繁但群间无联系 Q接近0（随机网络）：随机派对，每个人随机聊天，无法划分出明显的小圈子 Q为负值（反模块化）：刻意避免和“自己圈子”的人交流，反而只和“外人”互动在PRN中：高Q（如12Å截断）：蛋白质可清晰分成几个结构域（N叶、C叶），符合真实结构低Q（如16Å截断）：所有残基混在一起，失去模块边界，失去生物学意义重要说明：本文中使用modularity作为评估指标来量化网络的模块化程度，但论文并未详细说明具体的模块划分算法（如Louvain方法）或列出每个模块包含哪些残基。重点是通过观察modularity随截断距离的变化趋势（特别是在12-14Å范围内的斜率突变）来确定最优截断值，而不是深入分析模块的具体组成。边密度（Edge Density）实际边数与可能的最大边数之比，定义为： \(\rho = \frac{2|E|}{n(n-1)}\) 其中 $ E $ 是实际边数，$n$ 是节点数斜率变化的物理意义通过分析模块化和边密度随截断距离的变化曲线，发现12-14Å范围是最优的截断窗口：斜率变化标志着网络性质的转变：小截断（<10Å）：网络稀疏，模块化高但连通性差，斜率较陡（模块化随距离快速下降） 10-15Å范围：斜率明显变缓，这是从”模块主导”到”连通主导”的过渡区大截断（>15Å）：网络过度连通，模块化几乎消失，斜率趋平为什么斜率变化对应最优值：斜率最大处意味着网络性质变化最快，这是临界点在临界点之前：增加截断距离能够有效改善连通性，同时保持模块化在临界点之后：再增加截断距离只会模糊模块边界，不再带来新的结构信息二阶差分的数学意义：一阶导数 $f’(r)$：模块化随截断距离的变化率二阶导数 $f’‘(r)$：变化率的变化率（曲率）最大曲率点：一阶导数变化最剧烈的位置，即最优截断值插图显示：最大曲率出现在约12Å，因此确认其为最优值这一发现与先前LSP研究的经验一致，也符合蛋白质结构中邻近残基通常定义在12Å左右的常见做法。与传统LSP对齐方法的对比为了验证LSP-MD方法的可靠性，研究者将其与传统LSP对齐方法进行了系统对比：图6：LSP-MD与传统LSP对齐方法的结果对比。该图验证了LSP-MD方法与传统方法的一致性，同时展示了更高的计算效率： (A) 度中心性（DC）值的相关性：散点图对比了LSP-MD分析500帧和传统LSP分析100帧得到的DC值（均来自相同的10 ns PKA轨迹，5次重复）。横轴为传统LSP的DC值，纵轴为LSP-MD的DC值。关键发现：Pearson相关系数$r=0.91$，表明高度一致。大多数点沿对角线分布，误差棒（标准误差）较小，证明了LSP-MD能够重现传统方法的核心发现 (B) 介数中心性（BC）值的相关性：BC值的对比同样显示出显著相关性（$r=0.80$）。图中标注了三个具有高BC值的功能重要残基（K72、E91、D184），具体功能见下表 (C) 传统LSP的数据说明：图下方的说明文字指出，传统LSP方法由于计算复杂性限制，仅能分析轨迹的前100帧，而LSP-MD可以高效处理500帧。这种5倍的采样密度提升使LSP-MD能够更准确地捕捉热振动的统计特征使用相同的10 ns PKA轨迹，两种方法的效率和结果对比如下：对比维度 LSP-MD方法传统LSP对齐方法处理规模分析500帧仅能分析100帧（受限于计算成本）度中心性一致性 - $r=0.91$（与LSP-MD高度相关）介数中心性一致性 - $r=0.80$（与LSP-MD显著相关）关键功能残基的识别两种方法都识别出了一批具有高BC值的功能重要残基，具体如下：残基结构特征功能作用 K72 形成保守的K72-E91盐桥，连接N叶β折叠和αC-螺旋激酶活性态的标志，参与活性调控 E91 与K72形成盐桥稳定活性态构象，参与变构通信 D166 催化残基参与磷酸转移反应 D184 DFG基序的一部分参与镁离子结合和活性位点组织 F185 DFG基序的一部分其构象变化（DFG-in/out）是激酶活性的关键开关这些残基在PKA的功能和调节中发挥着核心作用，两种方法的同时验证确认了LSP-MD方法的准确性。方法的鲁棒性验证向量定义的独立性研究者测试了不同的残基向量定义对结果的影响（图S1），包括：标准向量：Cα→Cβ 长侧链向量：对于长侧链残基（如精氨酸、赖氨酸），使用Cα→侧链末端原子替代向量：对于甘氨酸，使用N→Cα或质心→Cα 结果显示，尽管不同向量定义导致绝对Z值有所差异，但DC和BC的相关系数均>0.95，证明中心性图谱对向量定义的选择不敏感。力场的独立性研究者使用两种不同的力场（ff14SB和CHARMM36）进行了对比模拟（图S2）。结果发现： DC相关系数：$r=0.98$ BC相关系数：$r=0.96$ 尽管两种力场对蛋白质动力学的描述存在差异，但LSP-MD捕获的中心性图谱高度一致，说明方法对不同力场具有鲁棒性。起始结构的独立性研究者从不同的起始构象（包括ATP结合态、抑制剂结合态等）开始模拟，并比较LSP-MD结果（图S4）。发现尽管局部动力学细节有所差异，但整体中心性图谱保持稳定，进一步确认了方法的可靠性。 Q&A Q1：LSP-MD方法与传统MD分析（如RMSD、RMSF、互相关分析）有什么本质区别？为什么要使用网络分析方法？ LSP-MD与传统MD分析方法的根本区别在于关注的物理量不同和信息抽象层次不同：表：传统MD分析方法与LSP-MD的对比 | 分析方法 | 关注的物理量 | 局限性 | 适用场景 | | — | — | — | — | | RMSD（均方根偏差） | 整体结构变化 | 无法区分局部稳定性差异，loop大运动和侧链小变化可能贡献相似的RMSD | 判断轨迹是否平衡、构象态聚类 | | RMSF（均方根涨落） | 单个残基涨落幅度 | 忽略残基间耦合关系，无法捕捉长程变构通信 | 识别高柔性区域 | | 互相关分析 | 残基间相关性 | 计算量大，相关系数矩阵难以直接转化为生物学洞察 | 初步识别残基间关联 | | LSP-MD | 残基对相对几何稳定性 | 需要构建PRN，计算复杂度略高于RMSF | 识别变构热点、分析局部刚性/柔性模块 | LSP-MD的独特优势聚焦相对几何：Z参数量化的是残基对的相对几何稳定性，而不是绝对位置变化。这对于识别局部刚性/柔性模块更为敏感网络化抽象：通过PRN将微观的几何涨落转化为宏观的中心性指标，天然地捕捉了多体耦合效应。高BC残基之所以重要，是因为它们位于多个通信路径的交汇处，这恰好对应了变构通信中的”热点” 物理意义明确：Z参数直接对应于构象熵（几何涨落越大，熵越大），而中心性指标则对应于该残基在变构通信中的重要性。这种从物理量到功能指标的映射链条清晰可解释一个形象的类比：想象一个城市交通系统。RMSD就像城市的整体繁荣度（所有人都在动），RMSF是每个人的忙碌程度（某些区域特别忙），互相关是人与人之间的联系矩阵（谁认识谁）。而LSP-MD的网络分析则识别出了”交通枢纽”——那些一旦堵塞就会导致全城瘫痪的关键节点。这些枢纽可能不是最忙的（RMSF不一定最高），也不是与所有人都有直接联系（度不一定最大），但它们位于不同区域之间的必经之路上（介数中心性高），因此对整体系统功能至关重要。 Q2：LSP-MD捕获的热振动（100 $\mathrm{cm^{-1}}$以下）与变构效应有什么因果关系？为什么这些微小振动能驱动远程的变构响应？这是一个深刻的物理生物学问题，涉及熵驱动变构的本质机制。Cooper和Dryden的理论预言可以通过LSP-MD方法得到直接验证，其物理逻辑如下：热振动的波数选择振动模式类型波数范围运动形式室温激发难易 LSP-MD覆盖高频模式 >200 $\mathrm{cm^{-1}}$ 键的拉伸和弯曲困难（能量高）否低波数模式 <200 $\mathrm{cm^{-1}}$ 扭动、剪切等集体运动容易（能量低）是 LSP-MD范围 5-100 $\mathrm{cm^{-1}}$ 侧链扭动、loop摆动充分激发完全覆盖熵-稳定性耦合机制一个残基对的热振动幅度（Z值）反映了其构象熵的大小。当配体在别处结合时，可能通过两种方式改变远程残基对的Z值：直接空间效应：配体的存在改变了局部空间位阻，远程残基的可动范围因此增大或减小间接网络效应：配体结合改变了某些关键残基（如铰链区残基）的稳定性，这种变化通过PRN传播，影响远程残基对的相对几何累积放大机制 Cooper和Dryden的理论框架提出，低波数振动模式（<200 $\mathrm{cm^{-1}}$）在生理温度下并未完全激发，可以在配体结合事件中被调制，从而导致构象熵的变化。单个残基对的熵变可能很小，但当多个残基对的熵变协同作用时，总效应可以被放大： \[\Delta S_\text{total} = \sum_i \Delta S_i\] 这种累积效应可以产生显著的自由能变化（$\Delta G = -T\Delta S$），足以驱动变构响应。许多变构调控的自由能差在5-20 kJ/mol范围内。从Z值到中心性的映射 LSP-MD的创新在于将微观的Z值通过PRN转化为宏观的中心性指标。高BC残基之所以重要，是因为它们连接了多个”熵变模块”。当这些模块的熵发生协同变化时，高BC残基就像是信息交换的枢纽，其稳定性变化会对整个网络产生放大效应。 Q3：10 ns模拟是否足以捕获所有与变构相关的热振动？对于那些发生毫秒级变构转变的蛋白质（如变构酶），LSP-MD方法是否仍然适用？这是一个关于时间尺度分离的重要问题，需要仔细区分不同类型的变构机制：时间尺度的层级结构蛋白质变构涉及多个时间尺度：时间尺度动力学过程捕获方法 LSP-MD应用皮秒-纳秒侧链热振动、loop快速摆动常规MD 直接分析微秒-毫秒构象态切换（open/closed）、domain运动增强采样MD 分态对比秒-分钟结合/解离、翻译后修饰生化实验/特殊方法不适用 10 ns的物理意义 LSP-MD聚焦于平衡态涨落，而非非平衡态转变。其假设是：蛋白质在特定功能态（如apo态或holo态）下，其热振动模式（由Z值分布表征）已经编码了该态的变构性质。如果两个态的热振动模式不同，那么其LSP-MD中心性图谱也应该不同。对于慢速变构系统的适用性对于那些发生毫秒级变构转变的蛋白质，LSP-MD的应用策略是：分别模拟不同功能态：对每个态（如open态和closed态）进行10 ns以上的模拟对比中心性图谱：计算两个态的DC和BC值，识别差异显著的残基识别变构热点：那些BC值在态间发生剧烈变化的残基就是变构通信的关键节点这种方法的物理基础是：即使构象转变本身很慢，但在每个态内部，热振动已经很快（皮秒-纳秒）达到了平衡。因此，10 ns模拟足以表征每个态的热振动特征，而态间差异则反映了变构效应。潜在局限与解决方案多亚稳态问题：如果10 ns轨迹在不同的亚稳态之间跳跃，Z值分布可能混合了多个态的特征。解决方案：使用聚类分析将轨迹分成不同亚稳态，分别分析构象异质性：某些蛋白质（如固有无序蛋白）本身就没有单一稳定构象。LSP-MD可能需要更长的模拟来捕获其系综特征。解决方案：使用多个短轨迹从不同起始构象开始模拟，构建综合的PRN 关键结论与批判性总结核心贡献物理基础明确：LSP-MD捕获的热振动波数范围（5-100 $\mathrm{cm^{-1}}$）与Cooper和Dryden理论预测的热激发模式阈值（<200 $\mathrm{cm^{-1}}$）高度吻合，为熵驱动变构提供了可量化的物理证据方法鲁棒性：中心性指标在不同模拟长度（图3）、采样率（图4）、向量定义（图S1）和力场（图S2）下保持稳定，证明方法捕获的是有意义的物理特征而非噪声截断距离优化：系统性地验证了12-14Å范围能产生最优的网络拓扑并保留关键结构信息（图5）与传统方法的连续性：LSP-MD保留了原始LSP对齐方法的核心结果（图6），同时用物理可解释的稳定性指标替代了ad hoc参数计算效率提升：这种连续性，结合改进的计算效率和更清晰的物理解释，使LSP-MD成为研究动力学驱动变构的实用可靠工具局限性与未来方向大尺度构象重排的挑战：一个悬而未决的问题是，LSP-MD记录的热动力学在涉及大尺度结构重排的变构系统中将如何表现。在这种情况下，局部熵特性可能在构象变化后发生改变。作者预期这些系统需要沿不同构象态分别取样分析。这些图谱的差异程度以及它们在什么时间尺度上达到平衡，仍有待确定。 BC的固有变异性：DC值的强相关性尤为重要，因为这一指标是研究熵驱动变构的主要关注点。相比之下，BC的相关性始终较低（图3F、4F、5D、6B、S1B、S2B、S4B）。这反映了BC的固有特性：它是依赖于最短路径的全局指标，边权重的微小变化就可能改变哪些残基被包含在这些路径中。因此，BC本质上比DC更易变，这是网络理论中公认的局限性。替代的中心性指标，如流介数（flow betweenness），可以应用于LSP衍生的PRN，但探索它们超出了这项以方法为重点的研究范围。小编锐评：本文基本上是在验证这个思想的可行性，各种指标什么的。显然不涉及大幅构象重排的变构过程，所以基本上是一个根据静态结构预测变构路径的增强版吧，可以作为未来工作流的一个步骤，比如边跑MD边根据这个工具修改CV？确实可能给DL训练提供数据？没验证是否适用于复合物，原则上应该可以吧如何对比两个体系，如ligand bound and unbound state，没给例子，似乎不好对比，只能各画各的图看不一样？
Molecular Dynamics · 2026-01-16
分子动力学揭示药物靶点变构通信路径：从动态网络到功能调控
 title: “MDPath：追踪蛋白质中的“悄悄话”——用分子动力学揭示药物靶点（如GPCRs）的变构通信路径” date: “2025-10-02” tags: [molecular-dynamics, sampling-and-analysis] — MDPath：追踪蛋白质中的“悄悄话”——用分子动力学揭示药物靶点（如GPCRs）的变构通信路径本文信息标题: MDPath：通过分子动力学模拟揭示药物靶点的变构通讯路径作者: Niklas Piet Doering, Marvin Taterra, Marcel Bermúdez, and Gerhard Wolber 发表时间: 2025年9月23日 (Accepted) 单位: 柏林自由大学生物、化学与药学系 (德国)，明斯特大学药物与医药化学研究所 (德国) 引用格式: Doering, N. P., Taterra, M., Bermúdez, M., & Wolber, G. MDPath: Unraveling Allosteric Communication Paths of Drug Targets through Molecular Dynamics Simulations. Journal of Chemical Information and Modeling. Published online September 23, 2025. https://doi.org/10.1021/acs.jcim.5c01590 源代码: https://github.com/wolberlab/mdpath 摘要理解蛋白质中的变构通讯对于基于结构的理性药物设计仍然是一个关键挑战。我们在此推出MDPath，一个用于分析分子动力学模拟中变构通讯路径的Python工具包，其核心是基于归一化互信息（NMI）的分析。我们以β₂-肾上腺素能受体、腺苷A₂A受体和μ-阿片受体为模型系统，展示了MDPath识别已知及新型GPCR变构机制的能力。该工具包揭示了β₂-肾上腺素能受体和MOR中配体特异性的变构效应，阐明了蛋白质-配体相互作用如何驱动构象变化。通过对ABL1激酶与变构和正构抑制剂复合物的分析，证明了该方法的广泛适用性。最终，MDPath为绘制蛋白质内部的变构通讯提供了一个开源框架，推动了基于结构的药物设计。背景变构（Allostery）是生物学中最基本的调控原则之一，它描述了一种“隔山打牛”的现象：蛋白质上一个位点的扰动（如配体结合或氨基酸突变）能够引起远处另一个功能位点的活性发生改变。这种远程调控使得药物分子不必直接作用于蛋白质的活性中心，而是可以通过结合在一个全新的“变构口袋”，来精细地调节蛋白质的功能，这为开发高选择性、低副作用的药物提供了巨大机遇。GPCRs、激酶等许多重要药物靶点都受到变构调控。然而，识别连接这两个远距离位点的“通讯线路”是一个巨大的挑战。这些线路并非静态的物理连接，而是由蛋白质内部残基间动态的、协同的运动所构成的复杂网络。静态的晶体结构往往无法揭示这些隐藏的动态信息，因此，分子动力学（MD）模拟成为捕捉蛋白质动态行为、研究变构机制不可或缺的工具。近年来，虽然涌现出多种用于分析MD轨迹以识别变构网络的计算工具，但它们大多关注于蛋白质整体的通讯网络，难以精确地分离出由特定配体结合所诱导的信号通路。此外，许多工具的设置复杂或并非开源，限制了其在药物研发领域的广泛应用。因此，亟需一个易于使用、开源且能系统性地、定量地描绘配体特异性变构路径的工具。关键科学问题如何从分子动力学模拟的海量数据中，系统性地、自动化地识别并可视化连接药物结合位点与功能远端位点的变构通讯路径？我们能否开发一个通用工具，不仅能确认已知的变构机制（如GPCR中的保守“微开关”），还能揭示配体特异性的调控网络（如激动剂和拮抗剂引发的不同信号通路），并为实验中观察到的突变效应提供合理的动力学解释？创新点发布MDPath开源工具包：提供了一个完整的、从MD轨迹分析到三维可视化的Python工具包，用于系统性地研究蛋白质变构通讯，其代码已在GitHub上开源。基于归一化互信息（NMI）的路径识别：采用NMI来量化残基间动态运动的相关性，并结合图论算法（Dijkstra）来寻找“信息流”最优的路径，为变构分析提供了数学上严谨且物理上直观的方法。配体特异性路径分析：实现了从特定配体接触残基出发追踪通讯路径的功能，能够清晰地区分不同配体（如激动剂与拮抗剂）引发的不同变构信号网络。广泛的验证与应用：在GPCRs和激酶这两大类重要药物靶点上成功验证了该方法，不仅重现了已知的保守变构基序，还为实验突变数据提供了新的机理见解。研究内容分子动力学模拟方法体系构建与参数化：研究使用了多个GPCR体系和ABL1激酶体系。GPCR结构来源于PDB数据库，包括激动剂结合态（β2：7DHI，A2A：2YDO，MOR：8EFQ）和拮抗剂/反向激动剂结合态（β2：5JQH，A2A：5MZP，MOR：7UL4），ABL1激酶结构为8SSN。所有体系使用MOE 2022.2进行预处理，包括缺失环区建模、突变回归野生型序列、添加缺失原子等。模拟软件与力场： GPCR体系：使用OpenMMDL进行体系构建，OpenMM进行MD模拟 ABL1体系：使用CHARMM GUI进行体系构建力场选择：蛋白质使用AMBER14SB力场，脂质使用Lipid21力场，配体使用GAFF2力场（ABL1体系中阿西米尼使用OpenFF）溶剂模型：TIP3P水模型，0.15 M NaCl离子浓度模拟参数：所有体系均进行能量最小化和0.5 ns平衡后，在NPT系综下运行3个独立的200 ns生产模拟。温度控制在300 K（Langevin动力学），压强控制在1.0 atm，时间步长2 fs，每个重复记录1000帧轨迹用于后续分析。核心方法论深度解析：MDPath的工作原理与流程 MDPath的核心思想是将蛋白质看作一个信息传递网络，利用MD模拟捕捉其动态行为，再通过信息论和图论的工具来寻找信息传递效率最高的“高速公路”。图5：MDPath用于变构通讯路径检测的主要工作流程。工作流程分为三个主要阶段：输入阶段接收MD模拟轨迹文件（PDB拓扑和DCD轨迹），可选择性添加配体相互作用位点等参数；分析阶段首先计算残基主链二面角运动，然后计算归一化互信息矩阵，接着构建网络图并使用Dijkstra算法寻找最大NMI路径，最后进行层次聚类识别核心通路；可视化阶段生成多种格式的输出文件，包括NGL view（Jupyter笔记本）、PyMOL脚本和STL文件（用于ChimeraX等软件）。 1. 相关性分析：从原子运动到信息网络为什么选择二面角？ MDPath选择监测每个残基的主链二面角（$\phi, \psi$）的动态变化，而不是Cα原子的笛卡尔坐标。这是一个关键的方法学选择。因为笛卡尔坐标会受到蛋白质在模拟盒子中整体平动和转动的影响，直接计算坐标相关性会引入大量虚假的、无物理意义的噪声。而二面角是内坐标，它只描述了肽链局部的扭转运动，与分子的整体运动无关。因此，基于二面角计算出的相关性更能反映蛋白质内部真实的构象变化和信息传递，信噪比更高。如何量化“通讯”强度？——互信息与NMI MDPath采用信息论中的归一化互信息（Normalized Mutual Information, NMI）来量化任意两个残基（X和Y）之间的“通讯强度”。首先，计算两个残基二面角运动之间的互信息（Mutual Information, MI）： \(MI(X,Y)=\sum_{x}\sum_{y}P(x,y)\log_{2}\left(\frac{P(x,y)}{P(x)\cdot P(y)}\right)\) 公式的通俗解释互信息衡量了知道一个变量后，另一个变量不确定性减少的程度，可以理解为两个变量之间非线性相关性的量度。 $P(x,y)$ 是联合概率分布，表示残基X处于状态x（某个二面角角度范围）且同时残基Y处于状态y的概率。 $P(x)$ 和 $P(y)$ 是边缘概率分布，分别表示X处于状态x和Y处于状态y的概率。如果X和Y的运动完全独立，那么 $P(x,y) = P(x) \cdot P(y)$，比值为1，$\log_2(1)=0$，MI为0。如果X和Y的运动高度相关，那么 $P(x,y)$ 会远大于 $P(x) \cdot P(y)$，比值大于1，$\log_2$项为正，MI值就高。然后，使用每个残基自身的熵（Entropy） $H(X)=-\sum P(x)\log_{2}(P(x))$ 对MI进行归一化，得到NMI： \(NMI(X, Y) = \frac{MI(X, Y)}{\sqrt{H(X)\cdot H(Y)}}\) 公式的通俗解释 NMI通过除以两个残基各自信息熵的几何平均值，消除了变量自身复杂性的影响。这使得NMI的取值范围被限定在0（完全无关）到1（完全相关）之间。一个高的NMI值意味着两个残基在动态运动上是高度协同的，即使它们在空间上相距很远，也表明它们之间存在一条有效的“通讯”通路。 2. 基于图的路径分析：寻找最优通讯路径计算出所有残基两两之间的NMI值后，MDPath将蛋白质抽象成一个网络图（Graph）。它将每个氨基酸残基视为一个节点（node），并在空间上邻近（< 5 Å）的残基之间创建边（edge）。关键的一步是如何利用Dijkstra算法。Dijkstra算法是图论中一个经典的最短路径算法，它寻找的是图中两点之间权重之和最小的路径。然而，我们的目标是寻找累积NMI值最大（即信息流最强）的路径。为了利用Dijkstra算法，MDPath进行了一个巧妙的转换：它将每条边的权重（weight）定义为与NMI值成反比的量（例如 $w = 1 - NMI$）。这样，NMI值越高（通讯越强），边的权重就越小。因此，在这个权重被“反转”的图中寻找“最短路径”，就等价于在原始概念中寻找“信息量最大的路径”。通过对所有可能的残基对运行该算法并筛选，MDPath便可描绘出蛋白质内部主要的变构通讯网络。 graph TD subgraph "输入阶段" direction LR A1["**MD模拟轨迹**<br/>拓扑文件PDB"] A2["**轨迹文件**<br/>DCD格式"] A3["**可选参数**<br/>配体相互作用位点<br/>分析参数设置"] end subgraph "相关性分析阶段" direction LR B1["计算所有残基<br/>主链二面角φψ轨迹"] B2["计算残基对间<br/>归一化互信息NMI矩阵"] B1 --> B2 end subgraph "路径分析阶段" direction LR C1["构建网络图<br/>残基为节点NMI为边权重"] C2["Dijkstra算法<br/>寻找最大NMI路径"] C3["层次聚类<br/>识别核心通路"] C1 --> C2 --> C3 end subgraph "可视化输出阶段" direction LR D1["**NGL view**<br/>Jupyter交互式"] D2["**PyMOL脚本**<br/>结构渲染"] D3["**STL文件**<br/>ChimeraX等软件"] end A1 --> B1 A2 --> B1 A3 --> B1 B2 --> C1 C3 --> D1 C3 --> D2 C3 --> D3 结果与分析 1. 模拟体系的质量控制：确保动力学轨迹的可靠性图S3-S5：激动剂结合的GPCR在200 ns模拟过程中的A100激活指数变化。 A100激活指数的计算原理：A100是一个专为A类GPCR设计的通用激活指数，基于五个关键的跨膜螺旋间距离计算得出。该指数通过机器学习方法训练，使用了大量微秒级分子动力学模拟数据和268个已发表的X射线晶体结构进行验证。A100指数的分类准确性在二态模型中达到94%（活性态）和99%（非活性态），在三态模型（包括中间态）中对活性态、中间态和非活性态的准确性分别为63%、81%和89%。在分析通讯路径之前，必须确保MD模拟本身是可靠的，即蛋白质在模拟过程中保持在预期的功能状态（活性态或非活性态）。作者使用A100激活指数来监测GPCR的构象状态（分数 > 0表示活性态，分数 < 0表示非活性态）。补充材料中的图S3-S5显示，在所有激动剂结合的体系中，A100分数在200 ns的模拟时长内基本都保持在0以上，表明模拟轨迹很好地维持了受体的活性构象，为后续的路径分析提供了可靠的数据基础。 2. 验证：识别GPCR中的保守变构“微开关” 图1：(A) 沙丁胺醇结合的活性态β₂-肾上腺素能受体的完整路径图。(B) 卡拉洛尔结合的非活性态β₂-肾上腺素能受体的完整路径图。(C) 热图显示了在所有三个模拟重复的前500条路径中，A类GPCR保守基序残基的参与情况。图中蓝色和紫色路径表示变构通讯路径，路径的粗细反映通讯强度。子图(D-H)详细展示了特定基序的路径：蓝色路径穿过CWxP基序(D)和PIF基序(E)，橙色残基标记关键基序位点。在非活性态中，蓝色路径通过NPxxY基序(F)和DRY基序的离子锁结构(G,H)。热图计算方法：图1C的热图统计了前500条最强通讯路径中每个保守基序残基的出现次数。对于每个基序（如CWxP、PIF、NPxxY、DRY），计算该基序内所有残基在路径中的参与频率，然后取该基序内任一残基的最大出现频率作为该基序的代表值。这种计算方式能够量化不同功能状态下各个保守”微开关”基序在变构通讯网络中的重要性。热图使用对数标度以更清晰地显示频率差异，颜色越深表示该基序在相应条件下的参与度越高。 GPCR的激活过程依赖于几个保守的氨基酸基序（”微开关”）的协同运动。MDPath的分析结果与已知的生物学机制高度吻合。在活性态受体(A)中，可以看到从细胞外域延伸到细胞内域的蓝色路径。非活性态受体(B)显示不同的路径模式。如图1C热图所示，在激动剂结合的活性态受体中，与激活相关的CWxP和PIF基序在通讯路径中的出现频率非常高。相反，在反向激动剂结合的非活性态受体中，与稳定非活性态相关的NPxxY和DRY基序则占据了主导地位。 3. 解释：为实验突变数据提供机理支撑图2：(A) 腺苷结合的腺苷A₂A受体中，从T88到W246的路径。(B) DAMGO结合的μ-阿片受体中，通过关键枢纽Y328的路径。图中蓝色路径表示变构通讯路径，橙色残基标记关键位点，黄色分子为配体。在A₂A受体(A)中，蓝色路径连接T88³·³⁶（橙色）到激活开关W246⁶·⁴⁸（橙色），展示从TM3到CWxP基序的直接变构通讯，解释了T88突变导致受体活性降低的机理。在μ-阿片受体(B)中，蓝色路径汇聚于关键枢纽残基Y328⁷·⁴³（橙色），该残基位于NPxxY基序上方，作为路径分布中心控制向细胞内结构域的信号传递。实验表明，在A₂A受体中将T88突变会显著降低受体活性。MDPath的分析（图2A）首次发现了一条从T88直达激活开关CWxP基序的变构路径，为该实验现象提供了清晰的机理解释。同样，对于μ-阿片受体（MOR），MDPath也发现Y328是一个关键的路径“枢纽”（hub）（图2B），与其实验功能的重要性相符。 4. 洞察：绘制配体特异性的通讯网络图3：β₂-肾上腺素能受体中的配体特异性路径。(A) 激动剂沙丁胺醇结合的活性态中的路径集群。(B) 反向激动剂卡拉洛尔结合的非活性态中的路径集群。图中展示了两种不同的变构路径集群：蓝色和红色路径代表两个主要的通讯集群，路径粗细反映通讯强度。黄色分子为配体（沙丁胺醇或卡拉洛尔），橙色残基标记参与路径的关键位点。在激动剂沙丁胺醇结合的活性态(A)中，路径主要汇聚到激活相关的PIF基序，显示出典型的激活信号传递模式。在反向激动剂卡拉洛尔结合的非活性态(B)中，路径模式完全不同，主要连接到稳定非活性态的NPxxY基序。值得注意的是，N312⁷·³⁹在两种状态下都不是主要路径的组成部分，表明其主要作用可能是配体结合而非功能调控。 5. 方法的稳健性与拓展应用模型完整性的重要性：补充材料中的一个关键负对照实验表明，如果人为地截断GPCR的一个重要胞内环（ICL3），MDPath分析出的路径就会变得模糊不清，甚至出现矛盾的信号（如在激活模拟中出现失活路径）。这证明了使用完整的、高质量的蛋白质模型进行MD模拟是获得可靠变构路径的前提。变构调节剂的影响：补充材料（图S7）还探究了钠离子和胆固醇等变构调节剂对通讯路径的影响。结果显示，这些调节剂的加入虽然会改变某些路径的权重（如增强了钠离子结合位点周围的信号），但核心的通讯通路模式保持不变，显示了变构网络的稳健性。在激酶靶点中的应用：图4：(A) ABL激酶与波舒替尼（紫色路径）和阿西米尼（蓝色路径）结合的完整视图。(B) DFG基序被变构路径稳定在DFG-out构象。(C) 远端T212残基作为正构路径的终点。图中紫色路径起始于正构ATP结合口袋（波舒替尼结合位点），蓝色路径起始于变构肉豆蔻酰口袋（阿西米尼结合位点）。两条路径都汇聚到自抑制性SH3结构域，但通过不同的机制。子图(B)显示蓝色变构路径如何稳定DFG基序（橙色）保持DFG-out构象，为阿西米尼的变构抑制机制提供分子基础。子图(C)展示远端T212残基（橙色）作为紫色正构路径的终点，解释了该位点突变如何影响ATP结合口袋抑制剂的活性。为了证明方法的普适性，作者将其应用于著名的ABL1激酶。MDPath成功识别出由正构抑制剂（波舒替尼）和变构抑制剂（阿西米尼）引发的两条截然不同的路径，并首次从动力学网络角度揭示了阿西米尼的变构抑制机制。 Q&A Q1: 这个工具对于药物研发的实际价值体现在哪里？ A1: MDPath的价值主要体现在以下几个方面：理解药物作用机制：通过可视化不同药物（如激动剂vs拮抗剂）引发的特异性通讯路径，可以深入理解其产生不同药理效应的分子基础。指导理性药物设计：识别出的路径上的关键“枢纽”残基，可以作为新的药物设计靶点，或者用于指导对现有分子的结构优化。解释耐药性突变：MDPath可以找到连接药物结合位点与远处突变位点的变构路径，从而解释为什么一个远端的突变会影响药物的疗效。发现新的变构口袋：通过分析整个蛋白的通讯网络，有可能识别出此前未被发现的、对蛋白功能至关重要的“热点”区域，这些区域可能成为全新的变构药物靶点。 Q2: MDPath的分析依赖于MD模拟，那么模拟的时长和质量对结果有什么影响？ A2: 这是一个非常关键的实际问题。模拟的时长决定了构象采样的充分性。本文使用了200 ns的模拟，这对于捕捉局部、快速的二面角运动是足够的，可以很好地分析处于一个稳定状态的通讯网络。但如果想要研究从非活性态到活性态的完整转变过程，这种慢过程就需要更长的模拟或结合增强采样方法。模拟的质量，如力场的准确性、体系构建的合理性，直接决定了轨迹的物理真实性。如果模拟本身不准确（如本文补充材料中ICL3截断的例子），那么从中分析出的任何“路径”都将是不可信的。因此，高质量、充分采样的MD模拟是MDPath分析成功的基石。 Q3: 论文中提到了对路径进行“层次聚类”，这一步的目的是什么？ A3: Dijkstra算法会找到成百上千条独立的“最优”路径。许多路径在空间上可能是高度重叠、非常相似的，它们实际上代表了同一条宏观的通讯“干道”。层次聚类的目的就是将这些相似的路径自动地分组归类。MDPath通过计算不同路径上残基原子坐标的距离来衡量路径的相似性，然后将相似的路径聚成一类。这样做的好处是，可以从纷繁复杂的数百条路径中，提炼出几条（如3-5条）最具代表性的、结构上不同的核心通讯通路（path clusters），如图3A中显示的红色和蓝色两条截然不同的路径。这极大地简化了结果的分析和可视化，让研究者能更容易地抓住主要的变构机制。关键结论与批判性总结核心结论本文成功开发并开源了一款名为MDPath的Python工具包，用于从MD模拟轨迹中系统性地识别、分析和可视化蛋白质的变构通讯路径。该方法以残基主链二面角的归一化互信息（NMI）为核心，结合图论算法，能够有效捕捉残基间的动态协同运动，并绘制出信息传递的最优路径。在GPCRs和ABL1激酶等多个重要药物靶点上的测试表明，MDPath不仅能准确识别已知的保守变构基序和激活机制，还能揭示配体特异性的信号通路。 MDPath的分析结果与实验突变数据高度吻合，能够为突变如何影响蛋白质功能提供合理的动力学机理的解释。潜在影响为药物研发领域的研究者提供了一个易于使用且功能强大的开源工具，有助于加深对药物作用机制的理解，并指导基于结构的理性药物设计。其“配体特异性”路径分析功能，为研究GPCR功能选择性、偏向性激动等前沿问题提供了新的计算视角。存在的局限性该方法目前仅考虑了主链二面角的信息，忽略了侧链运动和水分子等其他可能参与变构通讯的因素。路径识别的准确性依赖于MD模拟的充分采样。对于涉及大的构象变化的慢过程，可能需要更长的模拟或结合增强采样方法。路径的可视化和解读在一定程度上仍需要研究者的专业知识和判断。未来研究方向将侧链构象、水分子网络等更多维度的信息整合到NMI计算中，以构建更全面的通讯网络模型。将MDPath与马尔可夫状态模型（MSM）等方法结合，分析不同构象状态之间的转变路径。利用MDPath分析更大规模的MD数据库（如GPCRmd），进行高通量的变构机制探索。
Molecular Dynamics · 2025-10-08

Touch background to close